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Research  on  Deterministic  Methods  of 
Seismic  Source  Identification 

Report  No .  2  -  Summary 

The  objectives  of  the  research  conducted  were  to:  (1)  Develop  and 
perfect  methods  of  synthesis  of  theoretical  seismic  radiation  fields  at 
regional  and  teleseismic  distance  ranges  in  order  to  provide  a  basis  for 
developing,  testing  and  understanding  seismic  source  discrimination 
methods;  (2)  Develop  formal  inversion  procedures  to  be  used  to  obtain 
quantitative  estimates  of  physical  parameters  of  seismic  sources  and  the 
medium  and  to  use  these  estimates  to  determine  underground  explosion  yields, 
and  as  a  basis  tor  source  identification  and  discrimination;  and  to 
(3)  Develop  the  means  and  methods  to  analyze  seismic  data  for  source 
identification  and  propagational  medium  characteristics  for  event 
discrimination  and  yield  estimation  purposes. 

-■In  this  report  we  describe  methods  of  modal  synthesis  which  provide 
very  detailed  theoretical  seismograms  in  the  near  and  regional  distance 
ranges.  In  the  discussion  we  give  the  mathematical  basis  for  the  method, 
along  with  some  examples  of  its  use.  We  also  describe  a  sophisticated 
automatic  signal  analysis  program,  which  is  to  be  used  to  generate  discrimination 
variables  from  observed  data.  The  program  and  underlying  procedures  were  tested 
using  synthetic  seismic  signals  in  order  to  evaluate  its  performance  in  the 
presence  of  interfering  noise.  The  signal  analysis  program  was  also  applied  - 
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to  observed  seismic  data.  The  purpose  of  this  program  is  to  isolate 

particular,  important,  seismic  phases  (discrete  signals)  and  to  obtain 

accurate  estimates  of  the  spectral  characteristics  of  these  signals 

for  purposes  of  event  discrimination.  Based  on  its  performance,  on  a 

variety  of  observed  and  synthetic  seismic  data,  it  appears  that  the  method 

can  provide  much  enhanced  discrimination/event  identification  capability, 

over  currently  used  methods.  Finally  we  illustrate,  and  discuss,  the 

use  of  theoretical  event  magnitude  data  (m^  and  M^)  as  a  bases  for  event 

identification.  In  particular  we  show  that,  based  on  theoretical  m^ 

versus  Mg  relations  obtained  by  synthesizing  seismic  wave  fields  in  the 

teleseismic  distance  range,  we  can  obtain  good  estimates  of  seismic  source 

parameters,  in  particular  stress  changes  and  failure  zone  dimensions  for 

earthquakes  and  effective  stress  and  nonlinear  zone  dimensions  for 

explosions,  and  can  use  these  parameters  to  form  a  discrimination  available 

* 

(stress  change  divided  by  failure  zone  surface  area)  which  can  serve  as 
a  very  effective  means  of  discriminating  underground  explosions  from 
earthquakes.  As  a  supplement  to  this  report,  we  include  summary  discussions 
of  (1)  anomalous  radiation  from  explosions  and  (2)  the  effects  of 
attenuation  in  the  earth  on  seismic  signals. 


I.  Introduction 


The  objectives  of  the  research  being  conducted  are  to;  (1)  Develop 
methods  of  seismogram  synthesis  using  mode  superposition  and  related 
methods,  (2)  Finalize  the  theory  for  source  inversion  by  modal 
decomposition,  (3)  Determine  the  anelastic  characteristics  of  the  medium 
using  known  source  characteristics,  (4)  Interpret  seismic  event 
discrimination  in  terms  of  the  physical  properties  of  the  source, 

(5)  Establish  seismic  event  discrimination  methods  from  formal  inversion 
techniques,  and  (6)  Establish  regional  discrimination  techniques  based 
upon  physical  properties  of  the  source. 

In  this  report  we  emphasize  research  results  relating  to 
(1)  seismogram  synthesis,  (2)  automated  regional  and  teleseismic  signal 
analysis  for  discrimination  and  yield  estimation,  and  (3)  source 
property  inversion  and  discrimination  of  events  in  terms  of  inferred 
physical  characteristics  of  the  source.  In  addition  we  have  included, 
as  appendices,  discussions  of  the  status  of  our  knowledge  and  understanding 
of  anomalous  seismic  radiation  from  explosions  (Appendix  1)  and  a 
similar  review  of  regional  and  teleseismic  attenuation  effects  on  body  and 
surface  waves  (Appendix  2) .  Finally,  in  the  Appendix  3,  we  include  the 
complete  set  of  theoretical  ra^  versus  Mg  curves  generated  for  use  in 
source  property  inferences  (i.e.  average  stress  changes  and  rupture  or 
nonlinear  zone  dimensions).  These  latter  curves  form  the  basis  for  the 
formation  of  a  discrimination  approach  which  would  rely  explicitly  on 
differences  of  source  (physical)  properties,  as  inferred  from  m^  and  Mg 
magnitude  data. 
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II.  Modal  Synthesis  of  Seismic  Radiation  Fields 

Normal  Mode  Superposition  -  The  Theoretical  Basis 

In  the  following  sections  we  present  a  summary  of  the  steps  one  must 
take  to  compute  the  normal  mode  eigenvalues  and  eigenfunctions.  Since 
this  subject  has  been  treated  in  the  past  by  a  number  of  researchers,  we 
have  omitted  much  of  the  detailed  analyses  and  given  only  the  major 
steps.  Basically,  the  analysis  given  here  closely  follows  that  of  Ben- 
Menahem  and  Singh  (1972)  and  Harkrider  (1964)  who  show  how  to  syn¬ 
thesize  Rayleigh  waves  for  a  general  point  source.  We  do  present,  how¬ 
ever,  a  number  of  new  results  and  methods  that  constitute  the  basis  for  a 
successful  mode  representation  of  a  complete  synthetic  seismogram. 

Coordinate  Systems  and  Layer  Conventions 

In  order  to  easily  match  the  boundary  conditions  at  the  horizontal 
layer  interfaces,  we  use  a  cylindrical  coordinate  system  which  is  shown  in 
figure  1  along  with  the  layer  numbering  conventions.  Since  we  will  be 
writing  solutions  of  the  wave  equation  in  each  individual  layer  and  then 
matching  boundary  conditions  throughout  the  stack,  we  employ  both  a 
global  coordinate  system  and  a  set  of  local  coordinates,  each  relative  to 
an  individual  layer.  The  origin  of  the  global  coordinate  system  will  be  at 
the  free  surface  with  the  positive  z-axis  pointing  down.  The  origin  of  a 
local  coordinate  system  will  be  at  the  top  of  the  layer,  with  the  radial  and 
azimuthal  coordinates  being  the  same  for  all  of  the  coordinate  systems. 

We  distinguish  between  the  global  and  local  vertical  coordinates  by 
using  an  unsubscripted  or  unsuperscripted  z  for  the  global  coordinate 
and  a  superscripted  for  local  coordinates  where  p  is  the  layer  index. 
The  depth  of  the  bottom  of  the  p0*  layer  in  global  coordinates  is  z  -  h ^ 
whereas  the  thickness  of  the  pth  layer  is  (^\(z  =  h^-i)  +  We  use  this 

dual  representation  for  all  of  the  functions  of  z  as  well.  Whenever  a 
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function  of  z  appears  without  a  layer  superscript  it  is  understood  that 
the  argument  will  be  in  global  coordinates  and  whenever  a  layer  super¬ 
script  does  appear,  then  the  argument  of  the  function  will  be  in  local 
coordinates.  Thus  for  some  function,  /( z) 


/(z) 


\z 


=  h(p-D  +  ^(p) 


The  Basic  Solution 


We  follow  Ben-Menahem  and  Singh  (1972)  and  write  the  Fourier 
transformed  solutions  of  the  elastic  wave  equations  in  the  plh  layer  in 
terms  of  the  vector  cylindrical  harmonics,  P,  B,  and  C  as: 


f  u^(r,tf,{<*5;W,fc)  kdk  (la) 

m=— •»  0 

=  radial,  transverse  and  vertical  components  of  the  displace¬ 
ment  in  the  pth  layer 


TP(  (lb) 

m=— •  0  x 

=  radial,  transverse  and  vertical  components  of  the  traction 
across  a  horizontal  plane  in  the  p01  layer. 

Where 

‘U‘m(r.*,1?phu.k)  =  x^(((p);u,k)Pam(r,^\k)  +  (2a) 

+  hu.k  )Bam(r,iJ;k)  +  x^(^p);u,k)Cam(r,'dIk) 

and 

=  A'fei)(f<,>);w^)Pam(r.iJ:fc)  +  (2b) 

+  X&H&ho.k)  B^fr.tf*)  +  X&}{pho.k)C aro  (r,T?;fc). 

We  define  the  vector  cylindrical  harmonics  as  follows: 


0 

0 


Jm(kr)etm* 


[Pm(r,tf;fc)]  = 


(3a) 


[Bm(r.i?;fc )]  = 


dJm(k  r 

2-  mi) 

d(k  r) 

e 

r) 

deiTn* 

(k  r) 

d-d 

0 

r) 

deint* 

(k  r) 

d-d 

[Cm<r,i?;fc)]  = 


dJm[kr)  if 

"  “  d(k  r)  e 

0 


The  Greek  subscripts  are  used  here  as  tensor  indices  which  range  from 
one  to  three. 

We  wish  to  impose  the  boundary  conditions  of  continuity  of  displace¬ 
ments  and  tractions  at  layer  interfaces,  zero  tractions  across  the  free 
surface,  and  no  sources  at  z  =  Since  we  have  expressed  the  solution  in 
a  separable  form  and  since  the  vector  cylindrical  harmonics  are  orthogo¬ 
nal,  in  order  for  these  boundary  conditions  to  be  met  we  must  impose 
them  on  the  z  dependent  integrand  factors  in  equations  (2).  Thus  we  can 
express  the  boundary  conditions  by  the  following: 

(1)  Continuity  of  displacements  at  the  layer  interfaces: 

=  x^(C(p))  (4a) 

-W  fb>)=o 

for  1  <  p  ^  N 

(2)  Continuity  of  tractions  at  the  layer  interfaces: 


,>•(?) -0 

for  1  <  p  <  N 

(3)  Zero  tractions  at  the  free  surface: 

[xl'Ko)]  =  [oj 
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(4)  The:  Soinmerfeld  radiation  condition  (no  sources  at  infinity)  applies 
to: 

and  *&>({<">).  (4d) 

We  specify  a  point  source  in  layer  p  =  s  and  at  r  =  0  and  z  =  z0  , 
({•(*)  =  (5isl).  We  shall  solve  the  problem  of  forced  vibrations  by  adding  a 
particular  solution  to  the  homogeneous  or  unforced  solution  in  the 
source  layer.  The  form  of  the  particular  solution  will  be  the  same  as 
equations  (1)  with  the  only  difference  being  that  the  z  dependent 
integrand  factors  will  be  different  from  those  of  the  homogeneous  solu¬ 
tion.  Now,  the  particular  solution  exists  only  in  the  source  layer  and  in 
all  the  other  layers  only  the  homogeneous  solution  given  by  equations  (l) 
exists.  We  can  thus  express  the  total  solution  in  the  source  layer  as  fol¬ 
lows: 

(r.*.<<*>;«)  =  -??<*£  +  P^ni  )kdk  (5a) 

m=-«°  0 

and 

=  T7  (T™  +  pf^)kdk.  (5b) 

m =--•  0 

Where  the  particular  solution  integrand  factors  are  denoted  by  a  P  super¬ 
script  and  are  given  by: 

Pu£!(r  .O.&hu.k)  =  wtii^hu.kW^r.Kk)  +  (6a) 

+  '^^)(f{s);w.fc)BaTn(r.i?;*)  +  w$({lshu,k)Can(r,-d;k) 

and 

PT(£{r.#J*hu.k)  =  Wti(&';U'VPam<r.*;k)  +  (6b) 

+  ^)(^);«.fc)BOTn(r.17;fc)  +  W &J(<4* >:&>,*  )Cam(r,-d\k  ). 

The  boundary  conditions  given  by  equations  (4)  must  be  modified  at 


the  upper  and  lower  source  layer  interfaces  as  follows: 
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(1)  Continuity  of  displacements  at  the  upper  and  lower  source  layer 
interfaces: 

=  *&<0)  +  wfil(O)  (7a) 

*£n+,)(0)  =  *&(*<•>)  +  W&K&')  (7b) 

(2)  Continuity  of  tractions  at  the  upper  and  lower  source  layer  inter¬ 
faces: 

=  X£U0)  +  W&(0)  (7c) 

*&+,,(0)  =  *&(*<•>)  +  *&>($<•>)  (7d) 

We  now  define  two  sets  of  vertical  eigenfunctions  from  the  homogene¬ 
ous  solution  z  dependent  integrand  factors  corresponding  to  the  P-SV 
and  SH  waves  as: 


[“*(*)]  = 


x  lm(z  )/k 

xZm(z  )/fc 

Xln(*)/kz 
Xz n{z)/kz 


,  x3m(z)/k 

?*<*»  =  (8b) 
Similarly  we  separate  the  P-SV  and  SH  components  of  the  particular  solu¬ 
tion  z  dependent  integrand  factors  as  follows: 


[Bs£>(t(s)>]  = 


wfc K&')/k 

(s))/fc 

W&«{s))/kz 


p-w*)]  =  <9b> 

The  R  superscripts  correspond  to  the  Rayleigh  wave  solutions  and  the 
L  superscripts  correspond  to  the  love  wave  solutions.  In  this  paper  we 
will  only  consider  the  P-SV  case  and  will  develop  the  SH  case  in  a  separate 
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study.  With  this  in  mind  we  drop  the  R  superscript  in  order  to  minimize 
notational  clutter.  We  also  note  that  the  azimuthal  subscript,  m,  has 
been  dropped  from  the  eigenfunctions,  [£■],  since  they  are  m  indepen¬ 
dent. 


The  Eigenfunction  Solutions  and  the  Propagator  Matrices 

The  eigenfunction  solutions  for  a  single  layer  are  given  by  Ben- 
Menahein  and  Singh  and  are: 

[£-0>)«-<r>)]  =  (10) 

where  [A^]  is  a  4  component  vector  of  constant  coefficients  that  are 
adjusted  to  satisfy  the  boundary  conditions,  is  a  (4  X  4)  matrix 

which  is  a  function  only  of  the  horizontal  phase  velocity  and  the  layer 
elastic  parameters,  and  [C^]  is  another  (4  X  4)  matrix  which  is  z  depen¬ 
dent.  The  values  of  q^  and  q'^  are  the  P  and  S  wave  vertical  phase  velo¬ 
cities  and  are: 


n(p)  = 


Vl  -  (c/otfr*)2 


for  c  < 


(p) 


(11a) 


•(p)  - 


T7^rr\ for  c>a 

for  c  <  ffo) 


Vl  -  (c /&*))* 


for  c  >  (fo) 


(lib) 


i  V(c/^)z  -  1 

and  and  v'^  are  the  P  and  S  wave  vertical  wavenumbers  and  are  given 
by: 


„(*)  = 


u 

,(p) 


,-(p) 


V' 


u 

0>) ' 


The  C  matrix  can  be  expressed  as  follows: 


(12a) 

(12b) 
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where 


p{j>)(v{j>)  ,^(p)) 


pto)^) '{<*))  =  exp(^^5) 


Q<P)(U'(P)'{<P)) 


Q  IP >(„■&>), {•£*))  ] 
(13) 

(14a) 


Q(p)(l/*(p)<f0»))  =  exp^'k^). 

The  B  matrix  is  given  in  the  appendix. 


(14b) 


We  can  rewrite  equation  (10)  in  terms  of  initial  values  of  the  eigen¬ 
functions  at  the  tops  of  the  layers  and  in  so  doing  we  define  the  propaga¬ 
tor  matrix,  [A^]. 

[E<*)((<P))]  =  [Ab>><{<P>)][i?(*)(o)]  (15) 

where  from  equations  (10),  (13),  and  (14)  it  is  easy  to  show  that: 

[4(p)(£(p))]  =  [fi0»)(7?(j’),77'b>))]  [fi <»» )(7?(w ))]-*.  (16) 

The  propagator  matrix  allows  us  to  compute  eigenfunctions  at  the  bottom 
of  a  layer  in  terms  of  eigenfunctions  at  the  top  of  the  layer  so  that  in  glo¬ 
bal  coordinates: 

[£(hb>l)]  =  [A(*H*°,))][£(h(p-*>)J.  (17) 

We  can  start  at  any  layer  interface,  p,  and  propagate  the  eigenfunctions 
to  any  other  interface,  q,  (where  p  <  q)  by  applying  the  boundary  condi¬ 
tions  expressed  by  equations  (4)  and  by  repeating  equations  (17)  and  in 
so  doing  we  define  the  interlayer  or  total  propagator  matrix,  [A^]. 

0'(h(»>)]  =  [A^-p)]  [£(*<»))] 

The  interlayer  propagator  matrix  is  defined  as  follows: 


(10) 
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[A^'P)]  =  [^(9+P+l-0(^(9+P  +  l-0)] 

t~p*  t 


(18) 


The  Dispersion  Function,  Eigenvalue,  and  EUipticity  Computations 


In  order  to  compute  the  eigenvalues  we  need  to  apply  the  boundary 
conditions  given  by  equations  (4c)  and  (4d)  of  zero  tractions  across  the 
free  surface  and  the  Sommerfeld  radiation  condition  in  the  infinite  half¬ 
space  at  the  bottom  of  the  stack.  The  Sommerfeld  radiation  condition  in 
the  half-space  can  be  insured  by  setting  two  of  the  constant  coefficients 
in  [A^*]  to  zero,  so  that  the  upward  propagating  terms,  and  Q^N\  are 
eliminated  from  the  solution  of  the  ODE  in  the  half-space.  More  precisely, 
the  Sommerfeld  radiation  condition  requires  that: 


A}"’  =  AT  =  o 

and  thus, 


A"0 


=  [flW]-1  [£  (A  <*-»>)]. 


(20) 


By  using  the  interlayer  propagator  matrix  in  equation  (20)  and  extracting 
the  first  and  third  rows  we  can  define  the  (2  X  4)  D  matrix. 


0 

0 


=  [£<*>]  [P(/i<*>)] 


where 


(21) 


z>1<ri)  =  ([a(v)r% 
d$~1)  =  ([swr% 


when  p  =  N- 1  andj  =  1,2, 3, 4 


(22a) 


and 


t  ([awr‘)u  Ag 

i=l 

I 


D^  =  £  ([fi(W)]-')3*  Ay 

t*i 


ti 


when  N  - 1  >  p  SO  and  j  =  1,2,3, 4.  (22b) 


We  evaluate  equation  (21)  at  the  free  surface  and  apply  the  remain¬ 
ing  boundary  condition  of  vanishing  tractions  at  the  surface  which 
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requires  that: 

E3(h<~  °>)  =  £,4(h(0>)  =  0 
where  =  0  and  thus: 

D[?  D$  f£,(o)  fo 

Dzi  Dzz  Ua(0)  ]  =  [o  '  (23) 

The  Rayleigh  wave  dispersion  function,  RF(a>,fc),  is  equal  to  the  deter¬ 
minant  of  the  (2  X  2)  submatrix  of  [D]  shown  in  equation  (23)  and  clearly 
equation  (23)  can  only  hold  if: 

rF («.*)  =  D[?Djg  -  D Is Dtt  =  o.  (24) 

We  find  the  eigenvalues  by  varying  u  and  fc  until  the  dispersion  function  is 
zero.  For  a  flat  layered  medium  the  eigenvalues  will  form  discrete  curves 
in  the  (u.fc)  space,  which  are  the  Rayleigh  wave  dispersion  curves,  for 
values  of  c  less  than  the  S-wave  velocity  in' the  half-space  and  for  values 
of  c  greater  than  the  S-wave  velocity  of  the  half-space  the  eigenvalues  will 
form  a  continuous  spectrum.  We  consider  the  frequency  to  be  the 
independent  variable  and  the  corresponding  wavenumber  eigenvalues,  for 
a  particular  frequency,  w,  and  for  the  nih  mode,  to  be  Rfcn(o).  Thus, 

RF(u,Rfcn(cj»  =  o  (25) 

In  order  to  compute  the  ellipticity,  Rcn(w),  we  evaluate  equation  (23). 

,  .  _  <».B*n(»))  M»)) 

B  "  w,ftfcn(cj)  D  lZl(w.Rfcn(w))  ^22  (w,Rfc„(w))‘ 

The  Eigenvalue  Numerical  Instability 


So  far  we  have  followed  the  traditional  Thomson-Haskell  matrix  for¬ 
mulation  which  has  been  used  for  many  years  to  compute  Rayleigh 
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dispersion  curves  At  this  point  we  diverge  from  the  traditional  approach 
which  sutlers  from  numerical  instabilities  at  frequencies  consistent  with 
body  wave  synthesis.  Recently,  Abo-Zena  (1979)  has  described  a  method 
by  which  the  Rayleigh  dispersion  function  can  be  computed  without  loss 
of  precision  at  arbitrarily  high  frequencies  and,  briefly,  the  method  can 
be  described  as  follows. 

If  one  were  to  algebraically  expand  the  elements  of  the  D  matrix  in 
terms  of  the  transcendental  functions,  P  and  Q,  given  by  equations  (14) 
and  then  use  this  expansion  in  computing  the  dispersion  function  from 
equation  (24),  one  would  find  that  a  number  of  these  expansion  terms 
would  exactly  cancel.  Unless  this  cancellation  is  done  explicitly  before 
the  dispersion  function  computation  is  coded  in  a  digital  computer  pro¬ 
gram,  then  the  cancellation  will  be  done  numerically  in  the  chain  of  arith¬ 
metic  executed  by  the  computer.  Unfortunately,  these  terms  which  can¬ 
cel  grow  exponentially  with  frequency  and  at  high  frequencies  they  are  so 
large  that  the  terms  which  do  not  cancel  are  lost  in  the  computer  word 
roundoff.  This  is  the  source  of  the  numerical  instability  in  the  Thomson- 
Haskell  formulation  and  ...Abo-Zena  has  explicitly  canceled  these  terms 
so  that  they  never  appear  in  the  computer  code. 

In  order  to  do  this  in  a  straightforward  manner,  Abo-Zena  defines  the 
Y  matrix  which  we  give  here  in  terms  of  the  previously  defined  D  matrix. 

Yj?]  =  Df'DiS?  -  for  j.k  =  1,2,3, 4  (27) 

The  elements  of  the  (4X4)  Y  matrix  constitute  all  of  the  possible  sub- 
determinants  of  the  D  matrix  of  order  two  and  clearly  the  Y  matrix  is 
anti-symmetric.  We  also  note  that  from  equation  (24)  it  is  obvious  that: 

rF(«.*)  =  r£W)  (28) 


As  we  shall  see  the  Y  matrix  is  also  handy  for  computing  the  eigenfunc¬ 
tions.  The  Y  matrix  can  be  computed  recursively  at  each  layer  interface 
using  the  following  relation. 

[yW]  =  [y{p*»]  [Atoft***)]  (29) 

The  key  to  Abo-Zena’s  method  is  the  way  in  which  he  computes  this 
recursion  relation  and  we  refer  the  reader  to  his  paper  for  the  details. 

Although  Abo-Zena's  method  as  given  in  his  paper  works  very  well  to 
cure  the  dispersion  function  numerical  instability,  it  does  so  at  consider¬ 
able  expense  over  the  Thomson-Haskell  method  in  terms  of  computer 
time.  We  have  carried  out  the  algebra  given  by  Abo-Zena  further  and 
have  come  up  with  a  much  more  efficient  algorithm  which  also  retains 
the  numerical  stability.  We  omit  the  details  of  this  rather  tedious  deriva¬ 
tion  and  give  the  algorithm  in  the  appendix.  We  might  also  point  out  that 
although  Abo-Zena  claims  that  his  method  is  successful  when  using  single 
precision  arithmetic,  our  experience  in  using  his  method  on  a  broad 
variety  of  structures  and  frequency  bands  indicate  that  in  general  one 
would  be  advised  to  use  double  precision  arithmetic.  This  becomes  more 
important  as  the  structure  depth  or  the  number  of  layers  increase  and  it 
is  especially  important  when  computing  the  eigenfunctions. 

The  Eigenfunction  Numerical  Instability:  The  Problem  and  the  Cure 

In  order  to  synthesize  Rayleigh  waves  for  buried  sources  one  must 
compute  the  eigenfunctions  and  these  computations  also  exhibit  numeri¬ 
cal  stability  problems.  Normally,  one  would  find  an  eigenvalue  by  zeroing 
the  dispersion  function  as  computed  by  Abo-Zena's  method.  This  value  is 
then  used  to  compute  the  ellipticity  so  that  the  surface  values  of  the 
eigenfunctions  are  defined.  The  total  propagator  matrix  to  the  desired 
depth  is  then  computed  and  this  matrix  is  used  in  equation  (18)  to 
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compute  the  eigenfunctions  at  the  desired  depth. 

We  tried  this  straightforward  approach  and  found  that  in  many  cases 
the  eigenfunctions  became  numerically  unstable  and  the  problem  is  more 
pronounced  as  the  frequency  is  increased.  Basically,  what  is  happening  is 
that  small  roundoff  errors  in  the  four  eigenfunctions  are  effectively 
amplified  by  subsequent  propagator  matrix  multiplications  until  the 
errors  get  bigger  than  the  correct  values  of  the  eigenfunctions.  These 
errors  can  be  thought  of  as  unstable,  drifting  errors  with  respect  to 
depth  since  there  is  no  ipherent  mechanism  in  the  computations  to  sta¬ 
bilize  them.  Thus  the  computed  eigenfunctions  at  half-space  depths  gen¬ 
erally  no  longer  meet  the  Sommerfeld  radiation  condition. 

In  order  to  illustrate  this  problem  we  use  as  an  example  the  only 
structure  for  which  we  could  readily  obtain  an  exact  analytic  expression 
for  the  eigenfunctions,  namely  an  infinite,  homogeneous  half-space. 
Using  the  propagator  matrix  from  the  appendix  and  replacing  the  hyper¬ 
bolic  trigonometric  functions  with  the  exponential  functions  P  and  Q,  we 
can  write  the  first  displacement  eigenfunction  as  follows. 


E\(z)  =  (  -(7-1)  +  *£i<n)  +  zp{z  )  (  “(T-1)  ~  K^iPr)  + 


.  9(z)  /  (7-1)  ^  .  1 

Where  we  have  dropped  the  layer  index  and, 


<  7  + 


(30) 


F(z)  =  cxp(— -  <p) 
RC1 


Q(z)  =  exp(^-/) 

Rc  1 

and  h Et  's  the  ellipticity  for  the  fundamental  mode  and  rcj  is  the  eigen 
phase  velocity  which,  of  course,  is  frequency  independent.  Since  the 
value  of  jjcj  is  less  than  the  S-wave  velocity  of  the  half-space,  the 


functions  <p  and  <p'  arc  real  and  positive  and  thus  the  exponential  func¬ 
tions  P  and  Q  have  real  and  positive  arguments.  If  -we  replace  the  ellipti- 
city  -with  an  analytic  expression  in  <p,  <p\  and  y  we  find  that  we  can  reduce 
equation  (30)  to: 

(3D 

The  cause  of  the  eigenfunction  numerical  instability  can  be  seen 
when  one  compares  equation  (30)  with  equation  (31).  The  Sommerfeld 
radiation  condition  requires  that  the  growing  exponential  solutions  vanish 
which  of  course  is  the  case  with  equation  (31).  This  comes  about  because 
the  terms  which  multiply  the  exponential  functions  P  and  Q  in  equation 
(30)  are  identically  zero.  In  order  to  show  this  one  must  substitute  an 
explicit  analytic  expression  for  the  ellipticity  which  in  general  is  impossi¬ 
ble  to  derive  but  for  this  simple  case  is  easy  to  derive.  If,  however,  one 
were  to  code  equation  (30)  in  a  computer  program  to  compute  the  eigen¬ 
function,  then  the  terms  which  multiply  the  growing  exponentials  would 
be  computed  numerically  and  instead  of  being  identically  zero,  they 
would  be  of  the  order  of  the  computer  word  roundoff  error.  These  small 
but  finite  terms  would  then  be  multiplied  by  the  growing  exponential 
functions  and  at  some  depth  the  eigenfunction  error  terms  would  over¬ 
come  the  correct  decaying  solutions.  Since  the  arguments  of  the  func¬ 
tions  P  and  Q  are  directly  proportional  to  frequency,  this  numerical 
becomes  more  pronounced  for  a  given  depth  as  the  frequency  increases. 

One  can  sec  from  equation  (30)  that  the  type  of  algebraic  cancella¬ 
tion  used  by  Abo-Zena  to  eliminate  the  dispersion  function  instability  can 
only  be  realized  with  the  eigenfunction  computations  if  an  explicit  ana¬ 
lytic  solution  for  the  ellipticity  can  be  derived.  In  general  such  an 
expression  cannot  be  derived  and  so  in  general  wc  are  stuck  with  using 
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equation  (30)  or  a  more  complicated  version  for  a  buried  layer. 

We  were  able  to  solve  this  problem  by  adding  an  additional  constraint 
on  the  eigenfunction  computations.  From  equations  (21)  we  can  see  that 
they  constitute  a  linear  system  of  two  equations  in  the  four  eigenfunc¬ 
tions  and  thus  we  can  solve  for  two  of  the  eigenfunctions  in  terms  of  the 
other  two.  We  can  express  the  two  traction  eigenfunctions  in  terms  of  the 
two  displacement  eigenfunctions  as  follows: 


Ea(h^]) 


y(p)  yb>) 

1  U  1  24 


EAhto)) 


[E4(hW)  j  “  "  y$F  -Y&  )  U 8(h°,))  J  (  ; 

In  our  search  through  the  literature  it  appears  that  Duncan  (1965)  came 
up  with  a  similar  relation  to  equation  (32),  but  he  did  not  use  it  in  the 
manner  we  have  here  to  control  eigenfunction  stability.  This  relation  can 
be  used  to  compute  the  ellipticity  in  terms  of  the  elements  of  the  Y 


matrix  which  gives, 


Rtn(u)  =  ~ 


W?  («.!*.(«))  r  {?<«**»(«)) 


Knv '  YgMM)  <«**•(«)>'  v  7 

In  order  to  apply  equation  (32)  to  the  half-space  example  problem  we 
must  divide  the  half-space  into  a  number  of  pseudo-layers  all  of  which 
have  the  same  elastic  parameters.  The  procedure  we  follow  to  compute 
numerically  stable  eigenfunctions  is  as  follows: 

(1)  Wc  evaluate  the  Y  matrix  at  each  pseudo-layer  interface. 

(2)  The  ellipticity  is  computed  from  equation  (33)  which  defines  the 
surface  values  of  the  eigenfunctions. 

(3)  We  compute  the  two  displacement  eigenfunctions  at  the  next 
pseudo-layer  interface  using  the  first  two  rows  of  the  pseudo¬ 
layer  propagator  matrix. 
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(4)  Equation  (32)  is  used  to  compute  the  two  traction  eigenfunctions 
at  the  pseudo-layer  interface  instead  of  using  the  last  two  rows 
of  the  propagator  matrix. 

(5)  Steps  3  and  4  are  repeated  to  the  desired  depth. 

The  reason  this  method  works  is  that  implicit  in  equation  (32)  are  the 
interface  boundary  conditions  and  the  Sommerfeld  radiation  condition  at 
the  bottom  of  the  structure.  Assuming  a  small  error  exists  in  the  dis¬ 
placement  eigenfunctions  which  could  come  about  due  to  the 
amplification  of  roundoff  errors  by  the  growing  exponential  propagator 
terms,  then  equations  (32)  effectively  introduce  small  compensating 
terms  in  the  traction  eigenfunctions  so  that  the  factors  which  multiply 
the  growing  exponential  terms  stay  small.  Without  such  a  correction  the 
error  terms  grow  exponentially  with  depth  until  they  overpower  the 
correct  solutions.  This  method  introduces  a  "feedback"  mechanism  into 
the  computations  which  controls  the  error  and  insures  that  it  stays  small 
with  depth. 

We  graphically  show  how  well  this  method  works  in  figures  2  and  3. 
The  half-space  problem  we  have  discussed  is  shown  in  figures  2a,  2b,  and 
2c,  which  are  plots  of  the  vertical  displacement  and  normal  traction 
eigenfunctions  as  a  function  of  depth  at  a  frequency  of  five  hertz.  The  S- 
wave  velocity  for  this  half-space  was  four  km/sec  and  Poisson’s  ratio  was 
0.25.  For  each  of  the  eigenfunctions  there  are  two  traces,  the  left  hand 
trace  is  the  computed  eigenfunction  plotted  on  a  linear  scale  and  the 
right  hand  tract;  is  the  difference  between  the  computed  trace  and  the 
exact  solution  from  equation  (31)  plotted  on  a  logarithmic  scale.  Figure 
2a  shows  the  numerical  instability  due  to  the  growing  propagator  term 
which  causes  the  eigenfunctions  to  "blow  up"  at  about  five  km  depth.  As 


can  be  seen  from  the  error  plots  this  instability  is  due  to  amplification  of 
the  initial  error  of  about  10-’6  which  is  the  double  precision  roundoff 
error.  In  figure  2b  we  have  placed  a  pseudo-layer  interface  at  five  km 
depth  and  then  applied  the  eigenfunction  stabilization  method  described 
previously.  Although  the  eigenfunctions  do  not  "blow  up"  they  clearly  are 
in  error  around  five  km  depth  and  they  are  showing  signs  of  instability 
once  again  at  about  ten  km  depth.  In  this  case  we  have  allowed  the  error 
to  become  large  before  we  applied  the  correction  at  five  km,  but  even  in 
the  case  where  the  eigenfunction  computations  are  obviously  in  error  the 
use  of  equation  (32)  forces  stability.  From  this  figure  we  can  see  that  the 
pseudo-layer  thickness  must  be  less  than  five  km  and  in  figure  2c  we  have 
used  2.5  km  thick  pseudo-layers.  In  this  case  the  error  never  gets  large 
and  the  use  of  equation  (32)  at  2.5  km  intervals  keeps  the  error  small. 
These  figures  show  that  this  stabilization  method  is  only  accurate  if  the 
error  is  not  allowed  to  become  large,  i.e.  the  pseudo-layer  thickness  is 
below  some  critical  value  which  will  be  inversely  proportional  with  fre¬ 
quency. 

From  a  practical  standpoint  we  will  not  be  very  interested  in  comput¬ 
ing  eigenfunctions  for  a  half-space  and  so  we  show  computed  eigenfunc¬ 
tions  for  a  more  typical  structure  in  figures  3a,  3b,  and  3c.  This  particu¬ 
lar  structure  is  two  five  km  thick  layers  over  a  half-space  and  the  eigen¬ 
functions  correspond  to  a  high  order  mode  at  five  hertz  frequency.  In 
these  figures  we  have  omited  the  error  plots  since  there  is  no  exact  solu¬ 
tion  to  compare  against.  Once  again  we  show  the  computed  eigenfunc¬ 
tions  using  the  normal  approach  in  figure  3a  and  as  with  the  half-space 
case  they  "blow  up"  at  about  five  km.  In  figure  3b  we  have  applied  equa¬ 
tion  (32)  at  the  structural  layer  interfaces  at  five  and  ten  km,  but  as 


before  we  might  suspect  that  the  errors  have  become  too  large  for  accu¬ 
rate  computations  of  the  eigenfunctions  at  and  below  five  km  depth.  We 
add  two  pseudo-layer  interfaces  at  2.5  km  and  7.5  km  and  the  results  can 
be  seen  in  figure  3c.  If  we  continue  to  add  pseudo-layers  to  the  structure 
and  make  the  pseudo-layer  thicknesses  smaller  the  resulting  eigenfunc¬ 
tions  do  not  change  from  those  shown  in  figure  3c  and  thus  we  conclude 
that  these  computed  eigenfunctions  are  good  approximations  to  the 
exact  solutions. 

The  Forced  Vibration  Problem  -  Source  Excitation 

The  source  z  -dependent  integrand  functions,  given  by  equations  9a, 
have  the  same  form  as  the  eigenfunctions,  given  by  equations  10,  11,  12, 
13,  and  14,  everywhere  except  at  the  source  depth  where  a  discontinuity 
in  at  least  one  of  the  functions  is  required.  We  follow  Harkrider  (1964) 
and  define  the  source  jump  vector,  [E^],  as  follows. 

[£&>]  =  rs£)(<5+{5))]  -  [S£,(«~(*))]  (34) 

where 

6  ■*"(*!  =  (5<*>  +  e 
=  <5(s)  -  e 

<5{s)  is  the  vertical  source  location  measured  from  the  top  of  the  source 
layer  and  e  is  a  positive  infinitesimal.  Using  the  source  jump  vector  we 
can  relate  the  source  integrand  components  at  the  top  of  the  source 
layer  to  those  at  the  bottom  by: 

[S&W*)}  =  (d«*)(^))][54f'(0)]  +  [A<SW’>-«5<S>)][E&>)  (35) 

Applying  the  boundary  conditions  as  with  the  eigenfunctions  we  can 
solve  for  the  total  solution  z -dependent  integrand  functions,  [Fi,0)],  at  the 


surface. 
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Du  D\l]  /''tm(O) 

])^  l)f}  >’zm(0)  =  ~  rfl<<>)l  \Z(z0))  '  [Em*]  ('^>) 

where 

[Z(*.)  ]  =  (37) 

and  is  the  total  propagator  matrix  from  the  surface  to  the  source  depth. 
The  total  solution  functions,  F,  become  eigenfunctions  when  the  source 
jump  vector  is  zero.  IVe  can  solve  for  the  surface  values  of  the  F  vector  in 
terms  of  the  Y  matrix. 


F,m(0) 

7bw(0) 


TE  M^n 


where 


A/lm  _ 

y(°) 

~  1  12 

0 

y(0) 

1  23 

y(0) 

■*  24 

M  2m 

0 

y(0) 

“■‘12 

y(°) 

*™  ■*  13 

y<°) 
"I  14 

We  can  now  write  the  integral  equation  for  the  displacements  at  the 
surface  due  to  a  buried  point  source.  From  equations  1,  3,  and  38  and 
only  considering  the  Rayleigh  contribution,  this  gives: 


RTZv(r.iS,0;cj)  =  S  /  ““ r)  -  r))eim*fcJ 

m  0  R*  k  r 

RU8(r,iS,0:to)  =  2  /  —  Jm(k  r)  eim*k  dk 
m  0  r 

•O 

RUa(r.a?,0;co)  =  1)  /  -An(*  r)  etm*  k'11  dk 

m  0  Ji" 


where 


(40a) 

(40b) 

(■10e) 


rF  =  Fl?(u.*) 

The  way  by  which  one  evaluates  the  integrals  given  by  equations  10  is 
one  of  the  major  topics  in  theoretical  seismology  and  in  this  paper  we  will 
approximate  the  integrals  with  the  residue  contributions  due  to  the  Ray¬ 
leigh  poles  along  the  real  wavenumber  axis.  Normally  this  approximation 
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would  give  good  results  only  for  the  surface  waves,  but  we  will  show  how 
this  approximation  can  be  made  to  give  good  results  for  bod)'  waves  as 
well  with  a  slight  modification  of  the  structure.  Thus,  considering  only 
the  residue  contributions,  we  can  write  the  solutions  to  equations  40  as 
f,  Tows. 


Rti  i(r,i?,0;cj)  =  —7 ii  £  £ 


M 


t  f-i  I®) 


a  y  r)  - 


dk 


RU2(r,i9,0;cj)  =  -rri  £  £ 


,,(2).  ,im 


imi)  j 


3fc 


|*=Rfcn(u) 

(41a) 


(41b) 


fc=Rfc„(u) 


RU3(r,i3.0;w)  =  -Tri  £  £ 


dRl> 


3fc 


(41c) 


fc=Rfcn(w) 


The  cylindrical  Bessel  functions  become  cylindrical  Hankel  functions  of 
the  second  kind  due  to  the  extension  of  the  integration  path  to  along 
with  the  symmetry  properties  of  the  integrands  (see  ). 

We  can  further  reduce  equations  41  to  a  compact  form  which  dees 
not  involve  the  propagator  matrix  to  the  source,  [Z],  but  instead  depends 
on  the  eigenfunctions  evaluated  at  the  source  depth.  This  was  done  by 
Harkrider  (1964)  and  we  omit  the  details  here. 


(42) 


HUa(r,'<J,z;a)  =  -Tri  £  S  An(u)  [E™  («)]  [tfn(z<>  ;u)]  R'kamn(r,i?,z  ;u)  Rfc„(cj)2 


where 


R^amn  (r.T?,  z  ,’cj)  =  Fin(z-,u)Vamn(r,i9;u)  +  BZn{z:u)Bamn{r,-0\u) 


(43) 
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An  («) 


8fc 
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fc  =Bfcn(CJ) 


[*£»(«)]  = 


v*  (s) 
^4m 


sr*(s) 

^lm 


l^=R^n(«) 


and  P  and  B  are  modified  vector  cylindrical  harmonic 


eigenvalues  and  are  given  by. 


[Pmn(r.**)] 
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H%\k  r) 


;fc=Rfcn(w) 
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3(fcr) 

H„(k  r)  de^*_ 
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k  =Rfc, 


(44) 

(45) 

evaluated  at  the 

(46a) 

(46b) 
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The  Locked  Mode  Approximation 

So  far  we  have  described  how  numerically  stable  spectra  oT  Hat  lay¬ 
ered  earth  normal  modes  can  be  computed  [or  a  point  source  at  arbi¬ 
trary  depth  and  for  arbitrarily  high  frequencies.  At  this  point  we  address 
the  problem  of  how  one  can  use  the  residue  contributions  to  approximate 
the  complete  integral  solutions  given  by  equations  (40).  We  start  by  first 
briefly  describing  the  properties  of  the  integrand  surface  in  the  complex 
wavenumber  plane. 

Many  researchers  have  investigated  the  wavenumber  integrand  func¬ 
tions  of  equations  (40)  (e.g.  Ewing,  Jardeztky  and  Press  [19??],  Aki  and 
Richards  [i960]).  These  functions  are  four  valued  due  to  the  dual  valued 
square  root  functions  of  phase  velocity,  ^ N '  and  ip'^N\  in  the  bottom  half¬ 
space  (see  the  Appendix  ).  Thus  the  Riemann  surface  will  be  four  leaved 
and  branch  points  occur  where  or  are  zero.  We  define  branch 

cuts  as  the  loci  along  which  the  real  parts  oT  p  and  ip'  arc  zero.  Figure  4 
shows  the  branch  points,  branch  cuts  and  pole  locations  Tor  real  fre¬ 
quency  and  assuming  no  attenuation.  Wc  have  omitted  the  mirror  image 
in  the  second  quadrant  since  the  integrand  topography  in  this  region  is 
not  important  for  evaluating  equations  (40). 

Our  approach  to  solve;  the  wavenumber  integrals  is  to  deform  the 
contour  of  integration  so  that  it  encircles  the  poles  and  then  use  the  resi¬ 
due  contributions  to  approximate  the  solutions  and  we  show  the 
deformed  contour  as  f  iri  figure  4.  We  can  replace  the  integrals  from  0  to 
+  «>  with  integrals  from  — to  +  »  by  replacing  the  Bessel  functions  with 
Hankel  functions  of  the  second  kind  (lludsonf  1969]).  The  contribution  of 
the  integration  path  at  \k\  =  «•  vanishes  so  that  the  wavenumber 
integral  equations  (10)  arc  equal  fo  the  residue  contributions  given  by 
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equation  (42)  plus  a  branch  cut  integral  contribution  which  can  be  seen 
in  figure  4. 

We  could  eliminate  the  branch  cut  integral  contribution  by  continu¬ 
ing  the  integration  contour  onto  the  other  Riemann  sheets  in  which  case 
the  contour  would  encircle  the  lower  two  quadrants  four  times  before  it 
would  close.  Wc  would  find  additional  residue  contributions  on  the  other 
sheets  which  would  be  located  off  of  the  real  axis  and  these  are  the  so 
called  "leaking"  modes  (Gilbert  [19??],  Aki  and  Richards  [1980]).  If  we 
were  to  include  these  leaking  mode  residue  contributions  we  would  have  a 
complete  and  exact  modal  solution  of  equations  (40)  without  having  to 
evaluate  any  line  integrals.  Unfortunately  these  leaking  mode  poles  are 
very  difficult  to  locate  and  without  the  development  of  some  efficient  and 
reliable  searching  scheme  these  residue  contributions  are  practically 
impossible  to  compute. 

If  we  ignore  the  leaking  mode  contribution,  the  resulting  synthetic 
seismograms  will  be  missing  the  early  part  of  the  coda,  namely  the  P- 
wave  arrivals  arid  the  early  S-wave  arrivals.  This  is  because  the  cutoff 
phase  velocity  is  the  S-wave  velocity  of  the  bottom  half-space  and  the 
modes  that  propagate  at  higher  velocities  are  the  leaking  modes  which 
we  have  neglected.  Wc  could  move  the  bottom  of  the  structure  down  for 
some  fixed  range  and  for  a  realistic  earth  model  this  would  result  in 
increasing  the  S-wave  velocity  of  the  bottom  half-space.  As  the  cutoff 
phase  velocity  increases  the  branch  point  would  recede  towards  the  ori¬ 
gin  in  the  complex  wavenumber  plane  and  leaking  mode  poles  that  were 
near  the  original  branch  point  location  would  migrate  towards  the  real 
axis  and  pass  through  the  branch  point,  as  it  recedes  to  appear  on  the 
real  wavenumber  axis.  These  modes  would  thus  be  converted  from  leak- 
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ing  modes  to  trapped  or  looked  modes  as  the  S-wave  velocity  of  the  bot¬ 
tom  half-space  increases  and  the  resulting  synthetic  seismograms  would 
contain  earlier  arrivals. 

The  essence  of  the  locked  mode  approximation  is  to  place  an  unreal¬ 
istically  high  velocity  half-space  at  the  bottom  of  the  structure  and  we 
refer  to  this  high  velocity  half-space  as  the  cap  layer.  We  then  require 
that  this  cap  layer  be  placed  at  a  depth  such  that  the  earliest  P-wave 
reflection  from  the  top  of  the  cap  layer  arrive  after  the  seismic  coda  of 
interest.  As  the  S-wave  velocity  of  the  cap  layer  increases  to  infinity  the 
entire  seismic  coda,  including  the  earliest  P-wave  arrivals,  can  be  syn¬ 
thesized  using  only  the  locked  mode  residue  contributions.  In  the  limit¬ 
ing  case  we  can  impose  a  boundary  condition  of  total  reflectivity  at  the 
top  oT  the  cap  layer  which  is  equivalent  to  specifying  zero  displacements 
along  this  interface.  This  changes  the  layered  half-space  problem  into  a 
layered  plate  problem  and  it  causes  the  wavenumber  integrand  function 
to  be  single  valued.  Thus  for  the  layered  plate  problem  no  branch  cuts 
exist,  there  is  only  one  sheet  to  the  Riemann  surface,  a  finite  number  of 
poles  are  located  on  Lhe  real  wavenumber  axis,  and  an  infinite  number  of 
poles  are  located  on  the  negative  imaginary  axis.  (We  might  point  out 
that  the  locked  mode  approximation  is  basically  the  same  as  the  box  nor¬ 
malization  method  which  is  well  known  in  quantum  mechanics  ([]).) 

There  are  a  number  of  problems  with  the  locked  mode  approximation 
which  wc  will  discuss  in  detail  in  following  sections  where  we  will  show 
with  numerical  results  how  and  under  what  circumstances  these  prob¬ 
lems  can  be  overcome.  The  problems  fall  into  three  categories  which  we 


list  below. 
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1.  Spurious  reflections  of r  of  the  cap  layer 

2.  Time  wrap  around  of  the  synthetic  seismogram  in  the  time  domain 
due  to  sampling  in  the  frequency  domain 

3.  The  truncation  effects  due  to  cuting  off  the  mode  sum  to  some  finite 
number  of  modes  (phase  velocity  filtering) 

Computational  Procedures 

We  have  implemented  the  locked  mode  method  as  described  in  this 
paper  on  a  DEC  POP  11/70  minicomputer  using  FORTRAN  computer  pro¬ 
grams.  This  computer  has  an  available  core  size  of  64 K  bytes  (l  byte  =  B 
bits)  for  program  instructions  and  64 K  bytes  for  data.  A  high  speed, 
large  volume  disk  is  also  necessary  in  order  to  handle  the  rather  large 
intermediate  data  files  (these  files  arc  sometimes  in  excess  of  several 
megabytes  in  size).  Single  precision  floating  point  variables  are  32  bits 
long  and  have  a  numerical  precision  of  about  8  decimal  digits  and  double 
precision  variables  arc  64  bits  long  with  a  precision  of  16  decimal  digits. 

Because  of  the  limited  core  space  for  both  data  and  program  instruc¬ 
tions,  wc  divided  the  computations  into  five  separate  computer  programs 
which  communicate  via  four  intermediate  data  files.  A  brief  description 
of  each  program  and  the  data  necessary  to  run  the  program  is  given 
below. 

1.  The  first  program  searches  for  poles  over  some  specified  phase  velo¬ 
city  range  and  for  some  fixed  set  of  frequencies.  Normally  the  fre¬ 
quencies  are  equally  spaced  from  zero  to  some  upper  cutoff  fre¬ 
quency  so  that  the  resulting  spectra  can  be  transformed  to  the  time 
domain  using  a  fast  Fourier  transform  program  without  using  any 
sort  of  spectral  interpolation.  The  other  input  parameters  are  the 
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structural  parameters  which  consist  of  the  P-wave  velocity,  S-wave 
velocity,  density,  and  thickness  of  each  layer.  We  assume  no  attenua¬ 
tion  at  this  point.  The  output  from  this  program  is  a  data  file  which 
consists  of  the  eigen  phase  velocities  of  all  of  the  modes  within  the 
prescribed  phase  velocity  window  and  at  the  prescribed  frequencies. 
In  order  to  avoid  numerical  problems  we  found  it  necessary  to  do  all 
of  the  arithmetic  in  double  precision.  The  eigenvalue  search  consti¬ 
tutes  the  major  expense  in  computer  time  and  it  also  required  the 
greatest  amount  of  time  to  develop. 

2.  The  second  program  computes  the  amplitude  term  given  by  equation 
(44),  the  group  velocity,  and  the  imaginary  shift  of  the  eigenvalue 
due  to  structural  attenuation  for  each  mode-frequency.  The  input 
eigenvalues  are  read  from  the  data  file  generated  by  the  first  pro¬ 
gram  and  the  only  other  inputs  are  the  values  of  Q  for  each  layer  and 
at  each  frequency  (it  is  quite  easy  to  account  for  frequency  depen¬ 
dent  Q).  We  assume  that  the  shift  in  eigenvalues  due  to  Q  can  be 
approximated  using  first  older  perturbation  theory.  Thus  the  real 
part  of  each  eigenvalue  remains  constant  and  the  first  order  shift  is 
purely  imaginary  (Archambeau  and  Anderson  []).  We  compute 
dRF/dc  numerically  (wc  found  t  his  to  be  the  cheapest  method)  at  the 
original  real  eigenvalue  and  Tor  real  P  and  S  wave  velocities,  we  then 
compute  the  dispersion  function  at  the  real  eigenvalue  but  using  the 
complex  P  and  S  wave  velocities  for  each  layer,  and  finally  we  apply  a 
first  order  Taylor’s  expansion  of  the  dispersion  function  to  compute 
the  imaginary  shift  of  the  eigen  phase  velocity.  Once  again  all  oT  the 
arithmetic  for  this  program  is  done  in  double  precision. 
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3.  The  third  program  computes  a  set  of  eigenfunctions  for  each  mode- 
frequency.  The  input  parameters  are  the  depths  for  which  the  eigen¬ 
functions  are  to  be  evaluated.  This  is  the  first  point  in  the  computa¬ 
tional  procedure  where  the  source-receiver  geometry  is  constrained 
since  the  source  and  receiver  locations  must  be  at  depths 
corresponding  to  the  computed  eigenfunctions.  All  of  the  arithmetic 
is  done  in  double  precision,  but  the  variables  on  the  output  file  are 
stored  in  single  precision.  We  found  that  subsequent  computations 
could  be  done  in  single  precision  without  undue  adverse  effects. 

4.  The  fourth  program  computes  the  receiver  displacement  spectra  for 
a  given  set  of  receiver  locations  and  source  parameters  using  equa¬ 
tion  (42).  Input  parameters  are  the  range,  azimuth,  and  depth  for 
each  receiver,  source  depth,  and  the  source  jump  vector  for  each 
mode-frequency  and  for  each  value  of  m  ,  the  azimuthal  index.  All 
three  components  of  the  receiver  displacement  vector  are  computed 
for  each  receiver  and  the  spectra  are  written  on  an  output  file. 
Phase  velocity  filtering,  group  velocity  filtering,  and  mode  number 
filtering  are  done  at  this  point.  All  of  the  computations  for  this  pro¬ 
gram  are  done  in  single  precision. 

5.  The  fifth  program  allows  for  seismic  instrument  convolution,  spectral 
amplitude  filtering,  and  computes  the  resulting  time  domain  signal 
for  each  component  of  each  receiver.  A  fast  Fourier  transform  sub¬ 
routine  is  used  to  transform  the  filtered  spectrum  to  the  time 
domain.  Graphics  software  is  used  in  this  program  so  that  the  result¬ 
ing  synthetic  seismograms  can  be  displayed. 

The  first  three  programs  are  basically  source  independent  (with  the 

exception  of  the  source  depth)  and,  for  a  given  structure,  only  the  last 
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III.  Seismic  Signal  Detection  and  Analysis  for  Event  Discrimination 

A  sophisticated  form  of  signal  detection/ 
analysis  has  been  developed  (Archambeau,  et  al . , 

1981a)  and  is  now  incorporated  in  a  new  computer 
program.  We  generally  refer  to  the  pro¬ 
cedures  involved  in  this  time  series  analysis  approach  as 
multiple-pass  dispersion  filtering  (see  next  section  for 
details).  The  entire  approach  is  based  on  quasi -harncnic 
decomposition  of  the  time  series,  using  narrow  band  filtering 
techniques,  and  is  described  quantitatively  by  Archambeau, 
et  al . ,  1981b.  The  approach  has  great  flexibility  in  that, 
when  three  component  seismic  data  is  available,  both  polari¬ 
zation  and  wave  number  filtering  can  be  carried  out  jointly 
with  the  dispersion  filtering  to  achieve  a  high  degree  of 
selectivity  for  signal  detection  and  a  large  amount  of  signal 
information  for  all  the  signals  detected.  For  regional  data, 
where  a  large  number  of  seismic  phases  arrive  close  to  each 
other  in  time,  it  is  not  only  desirable,  but  quite  necessary, 
to  isolate  individual  phases  and  obtain  good  (uncontaminated) 
estimates  of  the  spectral  content  of  each,  their  individual 
arrival  times  and  sense  of  first  motion,  and,  ideally,  also 
their  polarization,  dispersion  and  wave  number  vector  orienta¬ 
tions  as  functions  of  frequency. 

The  new  computer  program  developed  has  the  capability 
of  performing  the  entire  complement  of  signal  detection/ 
analysis  procedures  mentioned.  In  the  present  framework,  we 
are  employing  it  to:  (1)  generate  the  signal  data  needed  to 
define  discrimination  variables  and  location  parameters  for 
regional  events,  and  (2)  to  investigate  the  robust  character 
of  candidate  discrimination  variables,  using  both  real  and 
synthetic  seismic  data,  as  functions  of  noise  levels,  event 
type  and  magnitude  and  structural  variations.  Ultimately, 
we  also  expect  to  use  this  program  to  automatically  generate 
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all  the  required  discrimination  variables  to  be  used  in  a 
multivariate  discrimination  procedure  for  both  regional  and 
teleseismic  events. 

Figures  1  through  3  provide  an  example  of  the  testing 
of  this  program  for  signal  isolation  and  spectral  analysis 
using  regional  synthetic  seismograms.  The  objective  here  is 
to  verify  that  the  filtering  correctly  isolates  regional  type 
signal  pulses  in  the  presence  of  noise  and  obtains  proper 
spectral  estimates  for  each  known  signal  pulse  in  the  wave- 
train.  Figure  1  shows  predicted  ground  motion  as  a  function 
of  distance  from  the  source,  in  the  near-regional  distance 
range.  The  source  was  chosen,  for  this  test,  to  be  an  explo¬ 
sion  at  great  depth  (10  km) ,  so  that  the  synthetic  seismograms 
would  be  relatively  simple.  The  method  used  to  generate  the 
synthetics  is  the  locked  mode  approximation  devised  by  Harvey, 
1981.  The  phases  arriving  coherently  across  this  distance 
range  can  quite  easily  be  identified,  so  that  at  a  given  dis¬ 
tance,  it  is  possible  to  identify  particular  reflected,  re¬ 
fracted  and  direct  phases.  (Phase  identifications  can  also 
be  checked  by  theoretical  ray  tracing  through  the  model  veloc¬ 
ity  structure  used.) 

Figure  2  shows  the  results  of  adding  seismic  noise  to 
the  synthetic  seismogram  at  50  km  from  the  source  (top  trace) 
and  then  processing  the  resulting  time  series  through  the 
dispersion  filter  procedure.  In  this  operation,  only  pulse¬ 
like  signals  with  little  or  no  dispersion  (i.e.,  body  waves) 
are  selected  as  signals  of  interest.  The  procedure  involves 
ultra-narrow  band  filtering  and  generation  of  envelope  and 
instantaneous  phase  functions  from  which  the  energy  arrival 
times  (group  arrival  times)  and  the  spectral  amplitude  and  phase 
at  these  times  are  measured  and  associated  with  the  center 
frequencies  of  the  narrow  band  filters.  The  desired  undis¬ 
persed  pulses  are  then  selected  by  a  pattern  search  in  the 
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frequency- time-amplitude  parameter  space,  with  those  pulses 
having  the  proper  dispersion  and  preselected  frequency  depen¬ 
dent  signal-to-noise  threshold  level  being  selected.  Pulses 
are  selected,  their  spectra  subtracted  from  the  origin  data 
and  the  resulting  "reduced"  origin  spectrum  is  inverted  to  the 
time  domain.  The  entire  filtering  and  search  process  is  then 
repeated.  This  procedure  is  continued  until  no  new  signal 
pulses  are  found.  (This  iteration  procedure  is  the  origin 
of  the  term  "multiple-pass  dispersion  filtering.")  Pulses 
obtained  in  this  way  may  overlap  in  time,  and  it  is  possible 
to  separate  phases  that  are  time  shifted  by  only  a  fraction 
of  their  apparent  period.  On  the  other  hand,  if  the  pulses 
are  significantly  smaller  than  the  uncertainty  or  resolution 
time  of  the  narrow  band  filter,  then  they  are  automatically 
combined  together  (spectrally  added)  at  the  end  of  the  itera¬ 
tion  process  to  form  a  single  pulse.  The  individual  pulse 
spectra  so  isolated  are  then  inverted  to  the  time  domain,  as 
shown  in  the  lower  six  traces  in  Figure  2  .  In  the  figure 
these  may  then  be  summed  to  give  a  single  time  series,  as  in 
the  second  from  the  top  trace  in  the  figure,  which  can  be 
compared  with  the  original.  As  can  be  seen  from  this  compari¬ 
son,  the  detection  and  isolation  of  the  pulses  by  the  filter¬ 
ing  operations  appears  to  be  remarkedly  accurate.  This  can 
also  be  verified  by  comparing  the  filtered  output  results 
with  the  original  synthetic  in  Figure  1.  Close  examination 
shows  that  every  significant  pulse  was  detected  and  properly 
isolated  with  no  false  alarms.  (One  very  weak  arrival  was, 
however,  missed  in  the  detection  procedure,  presumably  due 
to  its  very  low  signal-to-noise  ratio.) 

Figure  3  shows  the  spectral  estimates  obtained  for 
each  of  the  pulses  detected.  That  these  spectral  estimates 
are  very  reasonable  estimates  of  the  true  spectra  is  apparent 
from  the  fact  that  when  inverted  to  the  time  domain,  they  re¬ 
sult  in  pulses  that  very  closely  match  those  of  the  pure  syn¬ 
thetic,  without  noise. 
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In  addition  to  the  dispersion  and  spectral  estimates  of 
the  isolated  pulses,  it  is  also  possible  to  obtain  the  "times 
of  arrival"  of  the  phases.  These  are  indicated  by  the  dashed 
lines  in  Figure  2,  and  the  procedure  is  described  and  illus¬ 
trated  in  a  following  section. 

Figures  4  and  5  illustrate  the  application  of  the 
procedure  to  a  small  regionally  recorded  earthquake  from 
Southern  California.  Only  the  first  three  seconds  of  the 
time  series  is  shown  in  Figure  5  and  within  approximately 
a  one  second  time  interval  five  separate  overlapping  signal 
pulses  are  detected.  The  separate  pulse  waveforms  are  shown 
in  the  lower  five  traces,  which  have  been  inverted  to  the 
time  domain  from  the  spectral  estimates  shown  in  Figure  5. 

The  sum  of  these  pulses  (second  trace  in  Figure  4  is  re- 
markedly  similar  to  the  first  second  of  the  event  recording. 

The  timing  lines  shown  are  computed  in  a  manner  designed  to 
indicate  where  the  maximum  of  the  early  arriving  energy  occurs, 
rather  than  at  the  signal  pulse  onset  time. 

Figure  5  lists  the  arrival  times  of  the  isolated 
pulses  and  it  is  evident  that  they  arrive  within  a  few  tenths 
of  a  second  from  each  other. 

It  is  evident  from  these  examples  that  we  are  in  a 
position  to  rapidly  and  accurately  generate  signal  data, from 
which  event  discrimination  and  yield  or  magnitude  results  may 
be  obtained.  In  particular,  we  clearly  can  measure  frequency 
dependent  magnitudes  over  a  wide  band  of  frequencies  for  any 
and  all  isolated  phases.  We  can  also  automatically  determine 
arrival  times  and  first  motion  data  for  all  of  the  phases 
quite  reliably.  We  can,  therefore,  investigate  spectral  dis¬ 
criminants  for  P  ,  P  ,  S  ,  L  ,  R  ,  etc.  and,  as  well,  obtain 

n  g  n  g  g 

complexity  measures,  fault  plane  solutions  and  location  and 
source  depth  estimates  using  all  the  phase  information  in¬ 
cluding  S-P,  pP,  PP,  etc. 
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In  this  regard,  we  plan  to  form  candidate  discriminants, 
such  as  those  mentioned  earlier,  and  systematically  test  them 
against  both  real  and  synthetic  regional  seismic  data. 

In  our  applications  to  synthetic  data,  we  expect  to  be 
able  to  verify,  in  considerable  detail,  the  applicability  and 
limitations  of  the  various  discriminants  in  a  wide  variety  of 
situations  by  varying  model  structure,  source  type  and  noise 
conditions  in  a  systematic  manner.  The  realism  of  this  approach 
can  be  gauged  by  noting  how  well  complex  regional  seismograms 
can  be  predicted.  For  example.  Figure  6  shows  a  layered  model 
representing  the  crust-upper  mantle  in  the  region  just  east  of 
the  San  Andreas  fault  in  Southern  California.  Figure  7  shows 
the  computed  P-SV  modes  that  are  associated  with  this  struc¬ 
ture,  by  means  of  a  phase  velocity-frequency  plot.  Seismo¬ 
grams  generated  in  the  regional  distance  range  from  300  to 
450  km,  from  a  shallow  (1  km)  depth  explosion  source,  are 
shown  in  Figure  8.  It  is  evident  that  the  commonly  observed 
complex  features  of  the  wavetrain  are  present,  as  well  as 
the  simpler  Pn  phase.  These  seismograms  are  extremely  similar 
to  real  data  observed  in  the  Southern  California  region,  and 
it  is  clear  that  most,  if  not  all,  of  the  important  charac¬ 
teristics  of  regional  event  data  are  represented.  Thus,  we 
can  reasonably  expect  that  definition  of  a  robust  discriminant 
based  on  tests  with  synthetic  data  of  this  sort  will  also  apply 
to  real  data. 
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Figure  8.  Regional  synthetic  seismograms  generated  by  mode 
superposition  for  the  Southern  California  area. 
The  lines  are  theoretical  travel  times  for  Pn  and 
Pq  based  on  simple  ray  theory.  The  complex  Pg 
wavetrain  results  from  many  higher  mode  contribu¬ 
tions.  Similar  theoretical  experiments  show  that 
a  "P  coda"  is  generated  when  thin,  low  velocity 
sedimentary  layers  are  present  near  the  surface. 
These  synthetic  seismograms  are  being  used  to 
define  and  test  discrimination  variables  and 
methods. 
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Automated  Signal  Detection,  Timing  and  Spectral 

Estimates  Using  Advanced  QHD  Methods 

Inherently  critical  to  the  realization  of  an  automated 
discrimination  (and  yield-magnitude  estimation)  procedure  is 
the  creation  of  a  fast,  sophisticated  signal  detection/isola¬ 
tion  procedure  and,  coupled  with  it,  related  procedures  for 
the  accurate  analysis  of  the  signal  data  for  spectral,  wave¬ 
form,  and  timing  characteristics-  The  function  of  this  sys¬ 
tem  would  then  be  to  automatically  provide  the  basic  data  to 
be  used  to  define  discrimination  variables,  such  as  those 
described  in  the  earlier  section. 

In  this  section  we  will  describe  a  currently  operational, 
but  newly  developed,  automatic  system  that  is  now  being  used  to 
obtain  the  required  signal  and  event  detections  and  the  basic  signal 
date  required  for  multivarient  discrimination.  The  new  system  uses 

multiple-pass  dispersion  filtering  as  a  principal  method  for 
the  identification  and  isolation  of  signals,  but  also  generates 
polarization  and  wave  number  data  if  three-component  seismograms 
are  available,  and  can  use  polarization  and  wave  number  filter¬ 
ing  jointly  with  dispersion  filtering  to  identify  and  extract 
signals  from  the  time  series.  The  entire  approach  is  based  on 
quasi-harmonic  decomposition  (QHD),  coupled  with  amplitude 
pattern  recognition  in  the  time-frequency  plane  that  can  be 
created  from  the  narrow  band  filter  output  used  to  perform  the 
decomposition  of  the  time  series  into  a  set  of  quasi-harmonic 
components.  Specifically,  energy  arrival  times  are  determined 
from  the  maxima  in  envelope  functions  generated  from  the  narrow 
band  filter  set  used.  At  these  times,  the  envelope  amplitude 
is  proportional  to  the  Fourier  spectral  amplitude  at  the  center 
frequency  of  the  filter  for  the  signal  (or  noise)  energy  arriving, 
while  the  instantaneous  phase,  computed  from  the  narrow  band 
filter  output  at  that  time,  is  directly  related  to  the  Fourier 
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phase.  Thus,  for  a  large  set  of  filters  with  different  center 
frequencies,  one  can  obtain  a  sample  of  the  spectrum  of  the 
energy  arriving  within  the  time  series  window.  This  provides 
spectral  data  as  a  function  of  time  throughout  the  time  series. 
Further,  when  two  or  more  displacement  components  are  avail¬ 
able,  then  the  phase  difference  and  amplitude  ratios  between 
the  two  components  can  be  used  to  define  and  measure  the 
polarization  and  azimuth/emergence  angle  (wave  number  vector 
orientation)  associated  with  the  wave  field. 

With  such  information  it  is  possible  to  search  for 
patterns  in  the  group  arrival  time-frequency-amplitude  param¬ 
eter  space  for  energy  arriving  with  a  particular  dispersion 
characteristic  and  the  required  polarization,  wave  number 
vector  and  signal-to-noise  ratio.  For  optimal  estimates  of 
the  spectral  and  group  arrival  times,  in  the  sense  of  maximum 
spectral  accuracy  and  resolution  in  time,  it  is  necessary  to 
prefilter  the  time  series  with  a  matched  filter  having  the 
inverse  of  the  sought  for  amplitude  and  dispersion  character¬ 
istics,  so  that  in  the  resulting  time  series  the  signals  of 
interest  are  pulse-like;  that  is  they  have  rather  flat 
spectra  with  little  or  no  dispersion.  In  this  case,  all  the 
signals  of  interest,  including  surface  waves  having  initially 
strong  dispersion,  would  have  pulse-like  character  and  the 
pattern  search  in  the  arrival  time-frequency  plane  would 
always  involve  a  search  for  undispersed  energy  arrivals  over 
the  frequency  band  covered  bv  the  filter  set. 
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An  important  feature  of  the  analysis  program 
is  a  multiple-pass  or  iteration  procedure.  In  this  case, 
pulse  arrivals  are  selected  throughout  the  time  series  using 
the  procedure  indicated  above.  Then  the  pulse  spectra  esti¬ 
mated  from  this  "first  pass"  through  the  data  are  subtracted 
from  the  spectrum  of  the  entire  time  series.  This  reduced 
spectrum  is  then  inverted  back  to  the  time  domain  to  produce 
a  "reduced"  time  series,  wherein  the  previously  detected 
pulses  have  been  removed.  This  new  time  series  is  then 
reprocessed  to  generate  new  time  varying  spectral  data  and 
the  data  is  searched  for  additonal  pulses  satisfying  the 
fixed  signal  criteria.  This  process  can  then  be  repeated  until 
all  "signals"  have  been  detected. 

Since  there  are  normally  overlapping  arrivals  in  the 
time  series,  it  is  common  for  one  of  the  arrivals  to  over¬ 
whelm  the  other  or  to  interfere  with  it  in  such  a  way  as  to 
make  it  violate  the  signal  selection  criteria  on  the  first 
pass.  However,  if  the  inteferring  signal  has  been  selected 
and  then  removed  (entirely  or  partially) ,  it  is  usually  the- 
case  that  the  second  signal  can  be  detected  and  properly 
identified,  because  of  the  reduction  in  interference.  Further, 
it  is  not  possible  to  obtain  totally  precise  spectral  esti¬ 
mates  of  any  signal  pulse  on  a  single  pass  through  the  time 
series,  and  so  when  the  reduced  spectra  for  the  entire  time 
series  is  inverted  to  the  time  domain  there  is  always  a 
residual  part  of  the  detected  signal  pulse  remaining.  This 
residual  is  usually  detected  as  "another"  signal  on  the  next 
pass  through  the  data,  if  it  is  above  the  background  noise 
level.  Hence  it  is  possible  to  recover  more  of  the  signal 
by  repeating  the  selection  process  on  the  reduced  time  series, 
as  well  as  to  recover  neighboring  signals  that  were  missed  due 
to  interference  effects. 
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Thus,  the  reduction  procedure  is  repeated  over  and 
over  again  until  no  pulses  satisfying  the  (fixed)  signal 
criteria  are  found.  At  this  stage  the  process  terminates. 
However,  because  of  the  possibility  of  repeated  detections 
of  the  same  signal,  as  just  described,  it  is  necessary  to 
consider  detection  associations  designed  to  combine  together 
(add)  the  spectra  of  pulses  detected  within  the  same  small 
time  window.  The  size  of  the  window  within  which  all  detec¬ 
tions  are  considered  to  be  the  same  signal  is  taken  to  be 
proportional  to  the  resolution  time  for  the  narrow  band  filters 
used. 

Once  all  the  associations  are  made,  then  the  remaining 
distinct  seismic  pulse  spectra  can  be  synthesized  into  the 
time  domain  individually,  to  give  each  pulse  waveform  in  com¬ 
plete  isolation  from  the  noise  and  other  pulses  detected. 

The  superposition  of  the  individual  pulses  then  should  give 
an  "uncontaminated"  seismogram,  which  can  be  compared  to  the 
original;  in  quantitative  terms  by  a  cross  correlation. 

Because  the  pulses  can  be  completely  isolated  in  the 
time  domain,  it  is  also  possible  to  time  them,  that  is  to  ob¬ 
tain  the  "time  of  first  motion."  The  first  motion  time  can 
be  obtained  by  a  variety  of  methods,  but  the  most  .stable  pro¬ 
cedure  is  to  form  the  individual  pulse  envelope  functions  and 
use  the  time  of  the  first  inflection  in  the  envelope  as  the 
pulse  arrival  time.  Further,  the  sense  of  first  motion  can 
also  be  simply  obtained  once  the  first  motion  time  is  deter-  , 
mined,  by  measuring  the  algebraic  sign  of  the  displacement 
immediately  after  the  first  motion.  These  operations  are 
included  in  the  QHD  program  and  are  designed  to  provide  auto¬ 
matic  timing  information  for  location  purposes  and  fault  plane 
solution  data  for  network  discrimination. 

Figure  9  illustrates  the  multiple-pass  dispersion 
filtering  operation  for  a  regionally  recorded  small  earth¬ 
quake.  The  multiple-pass,  reduced  time  series  are  denoted 
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(1)  through  (6)  in  the  figure,  with  (1),  of  course,  being  the 
original  seismogram.  The  first  arrival  from  the  event  is  near 
the  17  second  time  point  in  record  (1).„  The  multiple  passes, 
in  this  example,  were  designed  to  isolate  the  earliest  arriving 
pulses  having  nondispersive  character  —  that  is  having  body 
wave  characteristics.  Thus,  the  multiple  processing  will  first 
find  signals  at  the  left  of  the  figure  and  progressively  work 
to  the  right  until  all  signals,  satisfying  the  particular 
criteria  used  here,  are  found. 

The  dotted  lines  in  the  figure  denote  first  motion  times 
for  each  signal  pulse  detected.  Clearly  two  "false  alarms"  are 
obtained  during  the  first  pass,  since  they  occur  before  what 
we  know  to  be  the  first  arrival  from  the  event.  (These  can  be 
eliminated,  of  course,  by  adjustment  of  the  detection  criteria 
parameters.)  On  the  other  hand,  the  first  arrival  pulses  are 
properly  detected,  along  with  several  clear  late  arrivals.  The 
reduced  seismogram  obtained  by  removing  the  detections  obtained 
in  (1)  is  shown  as  trace  (2) ,  along  with  the  additional  pulse 
detections.  Similarly,  trace  (3)  is  the  reduced  seismogram 
with  all  the  detections  of  the  first  and  second  passes  sub¬ 
tracted.  Again  the  detections  obtained  are  shown,  with  the 
additional  detections  appearing  later  in  the  seismogram.  The 
remaining  traces  are  similarly  reduced  seismograms  and  show 
progressively  later  pulse  detections  as  the  early  signals  are 
reduced  to  noise  level.  (Since  the  seismograms  are  individually 
normalized  to  unity  for  the  largest  amplitude,  for  plotting 
purposes,  it  appears  that  the  noise  level  increases  from  trace 

(2)  to  trace  (6).  This,  of  course,  really  just  means  that  the 
signal  level  is  being  continuously  lowered  by  the  reduction 
process  toward  the  noise  level.)  The  process  automatically 
terminated  after  the  sixth  iteration,  with  no  remaining  sig¬ 
nals  detected. 

Figure  10  shows  an  expanded  time  scale  representation 
of  the  final  results  for  the  event  shown  in  Figure  9.  The 


Automated  pulse  decomposition  of  an  observed  re¬ 
gional  distance  range  earthquake  (top  trace)  using 
QHD  based  dispersion  filtering.  The  second  trace 
is  a  superposition  of  the  single  pulses  obtained 
by  the  decomposition  method.  The  isolated  pulses 
are  shown  in  the  lower  six  traces.  Also  shown 
with  the  single  pulses  are  their  individual  enve¬ 
lope  functions,  which  are  used  to  define  the  pulse 
arrival  times.  The  arrival  time  is  taken  to  be  at 
the  time  of  the  (first)  inflection  point  in  the 
envelope  function.  The  variation  of  the  signal 
pulse  amplitude  immediately  after  this  time  is 
then  used  to  determine  the  "sense  of  first  motion" 
of  the  pulse. 
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Figure  11.  Spectral  estimates  of  isolated  signal  pulses  from 
the  California  earthquake  shown  in  Figure  10-  The 
points  correspond  to  spectral  data  obtained  from 
the  narrow  band  filters  employed.  The  lines  repre¬ 
sent  interpolation-extrapolation  of  the  data  such 
that  there  is  automatic  smoothing  and  extrapolation 
to  flat  spectra  near  zero  frequency  (far  field 
representation  of  the  pulses) . 
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first  six  signal  pulses  after  the  first  arrival  are  shown, 
along  with  the  original  seismogram  (top  trace)  and  the  super¬ 
position  of  the  isolated  pulses  obtained  (second  trace  from 
the  top) .  The  envelope  functions  for  the  individual  pulses 
are  also  shown,  along  with  the  first  motion  timing  lines  that 
are  automatically  determined  from  the  envelope  inflection 
points. 

Figure  11  shows  the  spectra  for  the  first  three  pulses 
shown  in  Figure  10.  it  should  be  noted  that  the  program  auto¬ 
matically  extrapolates  the  spectral  data  in  such  a  way  as  to 
give  a  far  field  representation  for  each  pulse,  and  this  is 
not  directly  comparable  to  the  original  seismogram,  which  is 
band  limited.  However,  comparison  of  the  top  two  traces  in 
Figure  10,  nevertheless,  shows  strong  similarity  between  the 
original  and  the  reconstructed  seismograms. 

As  was  implied  earlier,  and  is  demonstrated  by  this 
example,  the  procedure  is  capable  of  separating  very  closely 
time  spaced  pulses.  This  is  evident  from  the  arrival  time 
differences  listed  with  each  signal  spectrum  in  Figure  11, 
where  signal  arrivals  are  only  separated  by  0.4  and  0.14 
seconds. 
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IV.  Theoretical  Basis  for  Event  Discrimination  Using  vs  H 

In  order  to  provide  a  theoretical  basis  for  m^  vs  type  discrimination 
we  have  generated  synthetic  seismograms  at  teleseismic  distances  from  earth¬ 
quake  and  explosion  source  models  in  representative  earth  structure  models. 

The  appropriate  way  to  do  this  is  to  generate  synthetic  seismograms  in  the 
time  domain  and  to  directly  measure  magnitudes  from  the  predicted  signals 
at  a  large  number  of  azimuths  and  distances  and  to  obtain  averages  in  the 
same  manner  as  is  done  in  generating  observed  magnitude  data.  We  have  used 
stress  relaxation  source  models  for  the  earthquakes.  (e.g.  Archambeau,  1964, 
1968;  Archambeau  and  Minster  197$)  and  explosion  models  based  on  numerical 
calculations  (e.g.  Cherry  et.al.  1972)  and  imperically  derived  explosion 
source  models  based  on  near  field  observations  (Mueller  and  Murphy,  1971). 

Since  the  earthquake  data  to  be  analyzed  will  represent,  in  general, 
events  occurring  in  both  oceanic  and  continental  tectonic  environments,  it 
is  necessary  to  obtain  theoretical  results  for  earth  structures  representing 
both  tectonic  continental  and  oceanic  upper  mantle  characteristics.  We 
have  used  the  continental  structures  for  tectonic  areas  obtained  by  Archambeau, 
Flinn,  and  Lambert  (1969)  as  well  as  a  modified  version  of  the  anelastic 
Q  structure  obtained  in  this  same  study  for  the  body  waves.  The  surface 
wave  anelastic  structure  was  adopted  from  the  MM8  model  of  Anderson  and 
Archambeau  (1964).  The  oceanic  velocity  model  used  was  based  on  that  given 
by  Toksoz  and  Anderson  (1962),  with  the  Q  model  modified  from  the  continental 
model  in  a  manner  consistent  with  the  differences  between  island  arc  oceanic 
(see,  Barazanzi  et.  al.  1975)  and  tectonic  continental  mantle  velocity 
distributions.  That  is,  the  low  Q  zone  was  made  to  correspond  to  the  low 
velocity  zone  in  the  model. 
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In  addition  to  the  necessity  of  generating  synthetic  time  series  for 
different  mantle  structures,  it  is  also  clearly  necessary  to  vary  the 
source  depth  and  the  failure  zone  orientation  with  respect  to  the  free 
surface,  inasmuch  as  these  variations  will  have  important  effects  on  the 
magnitudes  observed.  Thus  three  "types"  of  events  were  considered:  45° 
thrust  earthquakes,  normal  dip  slip  earthquakes  and  (90°)  strike  slip  earth¬ 
quakes.  Theoretical  seismograms  were  then  generated  for  these  basic  event 
types  at  five  separate  depths  (10  km,  15  km,  25  km,  35  km,  and  45  km)  in 
either  an  oceanic  or  a  tectonic  continental  environment.  In  addition  the 
rupture  velocity  was  varied  with  depth  since  all  the  events  used  a  rupture 
rate  equal  to  .8  of  the  local  shear  velocity,  which  was  different  in  the 
two  structures  used  and  different  at  different  depths.  Finally,  of  course 
since  we  wish  to  generate  theoretical  results  relating  body  and  surface 
wave  magnitudes  to  stress  drop  and  rupture  zone  dimensions,  both  of  these 
source  variables  were  varied  for  the  different  events.  In  particular, 
event  length  dimensions  of  1  km,  2.5  km,  5  km,  and  10  km  were  used  with 
each  event  type,  as  well  as  event  stress  drops  of  .001  kb,  .01  kb,  .1  kb 
and  1  kb.  All  these  variations  of  structure,  event  type,  depth  and  failure 
parameters  resulted  in  the  generation  of  a  very  large  number  of  synthetic 
seismograms  of  course. 

The  theoretical  event  magnitudes  m^  and  Mg  were  obtained  from  the 
synthetic  seismograms  using  the  conventional,  and  rather  imprecisely 
defined,  procedure  for  the  measurement  and  computation  of  magnitudes. 
However,  only  the  .05  Hz  surface  wave  magnitude  was  obtained,  in  view  of 
the  variability  of  Mg  with  frequency,  and  so  comparison  of  the  theoretical 
results  to  observations  requires  the  use  of  Mg  data  observed  at  .05  Hz. 

Since  this  has  now  become  essentially  standard  observational  procedure,  this 
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presents  no  particular  difficulties.  Similarly,  only  the  vertical 

component  M  was  generat  d,  since  this  is  now  the  standard  measurement, 
s 

but  the  radial  component  Mg  value  can  be  easily  obtained  by  multiplication 
of  the  vertical  amplitudes  by  the  appropriate  ellipticity  factor. 

The  theoretical  magnitude  data  generated  by  this  rather  massive 
computational  effort  is  illustrated,  in  one  particular  form  in  Figure  (12) 
through  (14). 

These  figures  show  the  theoretical  body  and  surface  magnitudes  for 

the  three  basic  event  types  as  functions  of  depth,  with  event  rupture 

dimension  as  a  parameter,  and  for  the  two  basic  mantle  structures  used . 

The  stress  drop  for  the  events  was  at  a  fixed  100  bar  level.  These  three 

figures  together  constitute  the  entire  theoretical  m^  and  Mg  set  for  the 

events  having  a  fixed  100  bar  stress  drop.  For  events  with  larger  or  small 

stress  drops  the  magnitude  scale  is  simply  shifted  linearly.  That  is,  since 

source  spectral  amplitude  scales  linearly  with  stress  drop, then  for  1000  bar 

stress  drop  events  the  magnitude  scale  would  change  by  the  (base  ten)  log  of 

the  ratio  of  the  1000  to  100  bar  stress  drop  magnitudes,  or  by  one  magnitude 

unit.  Hence  for  1  kb  stress  drop  events,  the  magnitude  scale  can  simply  be 

increased  by  one  unit  and  the  proper  predicted  m,  and  M  values  would  be 

b  s 

obtained.  Hence  these  figures  contain  essentially  all  the  magnitude  information 
generated. 

The  important  characteristics  of  these  results,  especially  for  the 
applications  intended,  are  first  that  the  body  wave  magnitude  shows  little 
variation  with  event  hypocentral  depth  and  essentially  no  variation  between 
continental  and  oceanic  structures.  This  is,  of  course,  highly  desirable 
for  purposes  of  source  property  inversion  and  event  discrimination  based  on 
these  parameters.  Second,  as  expected,  the  Mg  results  shown  systematic 
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differences  in  continental  versus  oceanic  structure  and  a  rather  strong 
variation  with  depth  in  the  hypocentral  depth  range  0  to  20  km,  in  both 
structures.  However  the  Mg  variation  with  event  depth  in  the  range  from, 
somewhat  less  than,  20  to  50  km  is  mild  and  the  curves  are,  in  fact,  very 
near  flat  at  a  constant  Mg  level.  This  is  very  important  for  source 
parameter  inversion  using  Mg  data,  since  it  is  highly  desirable  that  the 
data  be  strongly  dependent  on  only  the  source  parameters  of  interest 
(stress  and  rupture  dimension)  and  essentially  independent  of  depth  in 
particular,  since  hypocentral  locations  are  usually  uncertain  to  at  least 
20  to  30%.  Thus,  we  see  that  because  of  the  rather  fortunate  circumstance 
that  neither  Mg  nor  m^  values  are  very  sensitive  to  hypocentral  depth  in 
the  range  from  about  15  to  50  km,  we  can  expect  to  largely  avoid  uncertainties 
in  stress  estimates  arising  from  source  depth  uncertainties  for  events  in 
this  depth  range.  Thus  only  for  events  in  the  0-20  km  depth  range  is  it 
likely  that  depth  is  a  critically  important  factor. 

The  appropriate  form  of  the  theoretical  magnitude  results  is  as  a  set 
of  theoretical  curves  in  a  plot  of  m^  versus  Mg  value.  A  set  of  continuous 
curves  of  this  sort  can  be  easily  generated  by  interpolating  the  (m^,  Mg) 
results  obtained,  using  fault  dimension  and  stress  drop  as  parameters  for 
the  curves.  The  set  of  Figure  (15)  through  (18)  give  examples  of  such 
curves  for  45°  thrust  and  normal  dip  slip  events  along  with  observed  data 
from  the  Alaskan-Aleutian  Arc  region.  The  Appendix  3  contains  the  entire 
set  of  theoretical  m^  vs  Mg  curves  for  both  oceanic  (Island  Arc)  and  Tectonic 
Continental  origins,  for  all  the  source  types  over  the  source  depth  range 
0-45  km.  The  lines  labeled  1  km  5  km,  etc.,  are  the  loci  of  events  in  this 
Mg  vs  m^  plane  having  constant  rupture  dimension,  but  variable  stress  drop. 


The  curves  labeled  .001  kb,  .01  kb,  etc.,  are  the  loci  of  events  in  the 
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Figure  15.  Surface  wave  magnitude  Ms  versus  body  wave  magnitude  mb,  showing  event 
data  in  the  depth  range  0-10  km  from  the  Aleutian-Alaskan  region.  The  curves 
are  theoretically  predicted  for  thrust  type  events  in  this  depth  range  and  are 
used  to  estimate  stress  drop  and  rupture  dimensions  for  the  observed  events. 
The  solid  circles  are  explosions. 
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BODY  WAVE  MAGNITUDE  -mb 

Figure  16.  Surface  wave  magnitude  Mg  versus  body  wave  magnitude  mb,  showing 
event  data  in  the  Alaskan-Aleutian  region  for  the  depth  range  10-20  km. 
The  theoretical  curves  are  for  a  45°  thrust  event  at  15  km  with  variable 
stress  and  dimension. 
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Figure  18.  mb-Ms  data  and  theoretical  curves  for  dip  slip  type  events  in  the  Alaskan- 
Aleutian  region.  Some  of  the  smaller  events  indicate  high  average  stress  drops 
near  1  kilobar  and  small  failure  zone  dimensions.  The  depths  of  these  events  are 
unknown  but  because  of  their  location  on  the  oceanic  side  of  the  trench  they  are 
interpreted  as  being  shallow. 
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magnitude  plane  with  constant  stress  drop  but  with  variable  fault  dimension. 

An  important  aspect  of  all  the  constant  stress  drop  curves  is  their  convergence 

in  the  region  of  low  magnitudes.  Thus  the  separation  of  the  "grid"  defined 

by  these  curves  is  not  constant  and,  particularly  for  small  magnitude,  high 

stress  drop  events,  the  convergence  is  such  that  estimates  of  stress  will,  have 

large  uncertainty  even  for  rather  small  uncertainties  in  the  observed  magnitude 

values.  This  is  particularly  the  case  for  uncertainties  in  values  for 

events  in  this  magnitude  range.  For  events  with  depths  between  20  and  50 

km  and  with  M  values  above  about  4.5  however,  the  stress  estimates  should 
s 

be  quite  good  average  values,  with  only  modest  uncertainty  due  to  magnitude 
data  and  depth  errors  and  bias  (e.g.,  typically  20  to  30%  uncertainties). 

In  addition  to  illustrating  the  nature  of  the  theoretically  predicted 
grid  of  versus  Mg  curves,  from  which  estimates  of  stress  and  failure 

dimensions  can  be  obtained.  Figure  (15)  through  (18)  show  magnitude  data 
for  the  Alaskan- Aleutian  region  (Seismic  Region  1) .  The  data  has  been 
separated  into  groupings  with  respect  to  event  type  and  depth  in  order  to 
remove  biasing  of  the  source  parameter  estimates  due  to  differences  in  the 
magnitudes  for  events  of  different  type  and  at  different  depths.  The 
distribution  of  the  data  in  the  m^-Mg  plane  is  typicaly  of  data  world-wide 
in  that  in  the  small  to  moderate  event  magnitude  range,  the  average  stress 
drops  range  from  about  1  bar  to  about  1  kbar,  with  the  mean  at  about  60  bars. 
The  range  of  event  maximum  dimension  is  typically  from  slightly  less  than 
1  km  up  to  around  30  km,  with  only  the  large  surface  wave  magnitude  events 
having  dimensions  larger  than  this.  The  mean  event  size  in  this  data  is 
about  10  km.  As  can  be  seen,  nearly  all  of  the  events  are  in  the  surface 
wave  magnitude  range  from  4.25  to  6.5  while  the  rn4  values  are  essentially 


all  below  6.0. 
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An  explanation  of  the  stress  drop  range  is  based  on  the  fact  that 
the  values  obtained  by  this  method  are  average  stress  drops  corresponding 
to  a  spatial  average  of  stress  changes  over  the  entire  failure  surface, 
or  the  average  taken  over  a  characteristic  dimension  of  the  order  of 
the  mean  of  the  longest  and  shortest  wave  lengths  of  the  seismic  data 
used.  In  this  case,  for  m^  and  Mg  data,  this  averaging  dimension  would 
be  of  the  order  of  50  km  for  events  with  fault  lengths  larger  than  50  km 
and  of  the  order  of  the  fault  length  itself  for  the  small  events.  We 
observe  from  the  data  distribution  in  Figure  (15)  through  (18),  that 
most  of  the  very  high  stress  drop  events  are  of  small  dimension  and  that 
most  of  the  very  low  stress  drop  events  are  associated  with  large  faults. 
This  strongly  suggests  that  stress  drop  averages  over  a  small  fault  will 
be  much  more  strongly  influenced  by  a  high  stress  at  the  hypocenter  and 
reflect  this  high  stress  in  the  average  stress  drop  inferred  from  the 
m^-M^  approach,  while  for  the  large  faults  the  average  is  over  essentially 
the  whole  failure  region  and  the  high  static  stress  drop  occurring  at 
the  hypocenter  does  not  influence  the  average  proportionately  as  much, 
so  a  low  average  is  obtained.  However  the  data  show  that  the  range  of 
stress  drops  associated  with  events  of  fixed  dimension  is  about  one 
order  of  magnitude  and  this  range  may  be  due  to  strength  variations  or 
to  the  fact  that  the  size  of  local  stress  concentrations  can  vary 
considerably. 

The  theoretical  ra^  vs  Mg  curves  for  explosions  in  various  types  of 
media  have  only  been  partially  computed.  However  their  form  is  suggested 
by  those  computed  by  Archambeau  et.  al.  1974  and  by  the  explosion  data 
points  shown  in  Figure  15.  In  particular,  explosions  in  any  type  of 
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event  media  are  predicted,  and  observed,  to  occur  to  the  right  of  the 
earthquake  lines  and  data  shown  in  Figure  (15)-(18)  .  (That  is,  they 
are  predicted  to  have  relatively  higher  and  lower  Mg  than  earthquakes.) 

In  terms  of  discrimination  using  source  parameters,  explosions  have 
parameters  characterized  by  small  source  dimension  (L)  and  high  (effective 
or  apparent)  stress  drop  (Acs')  compared  to  earthquakes.  Therefore  an 
appropriate  discriminant  based  on  measured  source  parameters  is  Aa/A, 
with  A  representing  failure  zone  surface  area.  For  such  a  variable, 
one  would  expect  the  earthquake  population  to  define  a  region  with 
relatively  low  values  of  Ao/A  compared  to  explosions.  Certainly  this 
is  not  surprising.  The  current  theoretical  results  do  however  provide 
the  means  of  forming  such  a  discrimination  variable  and  using  it  as  a 
standard  method  for  event  discrimination. 
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Tectonic  Generation  of  Anomalous  Radiation  from  Explosions 

by 

Charles  B.  Archambeau 


Introduction 

Our  current  understanding  of  "tectonically  generated" 
seismic  effects  is  that  whenever  the  medium  is  stressed  in 
any  manner  (overburden,  tectonic  stresses,  etc.)  an  explosion 
which  creates  a  cavity  and  shatter  zone  will  induce  an 
"anomalous"  motion  in  the  medium,  which  is  due  to  the  relaxa¬ 
tion  of  the  initial  stress  in  the  medium  around  the  fracture 
zone.  Anomalous  seismic  radiation  is  defined  to  be  that  part 
of  the  seismic  field  that  does  not  arise  solely  and  directly 
from  the  isotropically  occurring  conversion  of  the  explosive 
shock  wave  into  an  elastic  compressional  wave.  That  is,  any 
part  of  the  field  not  corresponding  to  a  pure  isotropic  com- 
p>'essional  source. 

Not  all  of  the  observed  anomalous  radiation  from  explo¬ 
sions  need  arise  from  stress  relaxation  effects,  since  aniso¬ 
tropy  and  other  local  inhomogeneities  (including  the  free 
surface  of  the  earth)  can  give  rise  to  similar  effects. 
Kisslinger  (1976)  and  Bache  (1976b)  have  reviewed  the  likely 
processes  for  production  of  the  observed  anomalous  radiation. 
It  is  doubtful  that  the  seismic  perturbations  due  to  effects 
other  than  "tectonic"  are  nearly  sufficient  to  account  for 
the  observed  anomalous  radiation  (see  Toksoz  and  Kehrer,  197  2; 
Archambeau,  1972;  Bache,  1976a);  especially  the  anomalous 
radiation  associated  with  long  period  surface  waves. 

Physical  Processes  of  Stress  Relaxation 

The  explosion  induced  stress  relaxation  has  generally 
been  associated  with  tectonic  stress  relaxation,  but  it  should 
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be  emphasized  again  that  any  stress  in  the  solid  medium,  what¬ 
ever  its  origin,  will  cause  a  radiation  effect  (e.g.,  litho- 
static  stress) .  It  is  generally  agreed,  however,  that  stresses 
of  tectonic  origin  are  most  likely  to  be  responsible  for  the 
larger  effects  observed.  What  is  currently  less  clear  is 
whether  the  stress  relaxation  is  due  to  the  creation  of  a 
nearly  spherical  shatter  zone  (e.g.,  Archambeau,  1972; 
Archambeau  and  Sammis,  1970)  ,  or  whether  it  is  largely  due  to 
"triggering"  of  an  earthquake  —  that  is  associated  with  induced 
faulting  corresponding  to  material  failure  having  a  strongly 
asymmetrical,  or  linear,  pattern  (e.g.,  Aki  and  Tsai,  1972). 

In  the  latter  case,  shatter  zone  induced  relaxation  occurs, 
but  in  addition  failure  along  a  preexisting,  or  newly  created, 
long  linear  fracture  is  also  thought  to  occur  and  would 
generate  additional  "anomalous"  radiation,  especially  at 
the  lower  frequencies.  The  distinction  between  these 
mechanisms  has  importance  for  discrimination  as  well  as  yield 
estimation  in  that,  if  stress  relaxation  due  to  triggered 
faulting  occurs,  then  much  larger  perturbation  of  the  low  fre¬ 
quency  radiation  would  be  expected  than  would  be  the  case  for 
shatter  zone  induced  radiation.  For  triggering  then,  one 
might  expect  perturbations  that  would  make  accurate  yield 
determinations  using  Mg  highly  uncertain  —  for  example,  the 
surface  waves  from  triggering  of  a  modest  sized  thrust  earth¬ 
quake  could  easily  completely  cancel  or  overwhelm  the  explosion 
generated  surfaces  at  all  azimuths  in  the  period  range  around 
20  seconds.  Further,  the  Mg  value  for  the  explosion  plus 
earthquake  could,  in  many  cases,  be  earthquake-like,  making 
discrimination  by  m^-Mg  problematic.  The  size  and  nature  of 
the  effects  of  triggering  would  be  difficult  to  accurately 
predict  or  to  correct  for,  since  not  only  would  information 
regarding  stress  drops  be  necessary,  but  knowledge  of  the  loca¬ 
tion  and  orientation  of  the  fault  plane  and  rupture  rate  would 
be  required.  On  the  other  hand,  if  the  physical  mechanism  is 


I 
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simply  stress  relaxation  around  the  roughly  spherical  shatter 
zone  created  by  the  explosion,  then  it  is  likely  that  the  per¬ 
turbations  in  the  seismic  radiation  would  be  relatively  smaller 
at  low  frequencies  and  more  easily  predicted  and  corrected  for; 
in  that  the  origin  and  geometry  of  the  failure  zone,  as  well  as 
its  formation  rate,  would  be  known  quite  well.  Further,  the 
size  of  the  shatter  zone  is  predictable  and  related  to  the 
explosion  yield.  (Prediction  of  the  size,  or  radius,  of  the 
shatter  zone  is  not  normally  available  from  either  scaling  law 
or  empirical  data,  or  from  numerical  code  calculations.  How¬ 
ever,  it  could  be  obtained  in  these  ways.  In  this  regard, 
see  the  recommended  research  section  below.)  Finally,  the 
tectonic  and/or  lithostatic  stress  levels  in  the  shallow 
depth  range  for  average  explosions  can  be  reasonably  esti¬ 
mated,  especially  if  earthquakes  from  the  area  can  be  studied, 
and  a  correction  for  stress  relaxation  effects  on  surface 
wave  observations  could  then  be  made,  probably  with  reason¬ 
able  confidence.  (Such  a  correction  would,  of  course,  entail 
uncertainties  and  surface  wave  based  yield  estimates  would 
have  to  be  compared  with  body  wave  yield  estimates.)  In  any 
case,  research  would  be  required  in  order  to  develop  the 
basis  for  any  such  correction  and  to  determine  its  accuracy. 

Of  the  two  possible  mechanisms  for  tectonic  release, 
it  seems  most  likely  that  both  processes  have  occurred  in  the 
past.  Fowever,  the  simple  spherical  shatter  zone  relaxation 
process  has  probably  been  the  mechanism  giving  the  anomalous 
seismic  radiation  observed  from  most  underground  explosions, 
with  triggering  occurring  for  only  a  few  events  (about  five 
percent  or  less  of  the  total) .  This  estimate  is  based  on  a 
variety  of  physical  evidence  (i.e.,  aftershock  locations 
versus  time,  location  of  the  "anomalous  event",  near  field 
strain  observations,  etc.),  plus  the  fact  that  the  anomalous 
surface  wave  radiation  can  be  adequately  explained  by  spheri¬ 
cal  shatter  zone  induced  stress  relaxation  for  explosions  with 
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F  factors*  less  than  around  one.  (When  F  *  1,  the  seismic 
energy  of  tectonic  origin  is  about  the  same  as  that  from  the 
pure  explosion.)  Even  for  explosions  where  F  factors 
considerably  larger  than  one  are  observed,  it  is  not  certain 
that  large  scale  faulting  is  required.  In  any  case,  it  is 
likely  that  triggering  is  quite  rare. 


Characteristics  of  Observed  Anomalous  Body  and  Surface  Wave 
Radiation  from  Explosions  and  Predictions  from  Explosion 
Induced  Tectonic  Release  "  “ 

Some  typical  examples  of  the  nature  of  the  anomalous 
radiation  from  explosions  are  shown  in  the  following  figures. 
The  effects  are  most  pronounced  for  surface  waves  and  in  this 
discussion  such  effects  are  emphasized.  Nevertheless,  body 
waves  are  perturbed  as  well,  but  the  effect  on  body  wave 
magnitudes  m^  are  much  less  than  the  effect  on  the  surface 
wave  magnitude  ;  so  long  as  m^  is  measured  from  the  first 
cycle  of  the  P  wave  motion  or  is  measured  spectrally  by 
narrow  band  filtering  using  the  group  arrivals  within  1  to 
1.5  seconds  of  the  direct  P  wave  first  motion. 

An  example  of  the  effects  of  tectonic  release  on  the 
first  arriving  compressional  waves  from  an  explosion  is  shown 
in  Figure  1  (Archambeau,  et  al. ,  1974).  The  event  modeled 
is  Handley  (1.1  megaton)  at  a  teleseismic  distance.  The 
tectonic  release  mechanism  modeled  is  stress  relaxation 
around  the  roughly  spherical  explosion  generated  shatter  zone, 
where  the  initial  prestress  is  taken  to  be  65  bars  and 
homogeneous.  This  is  a  modest  prestress  level,  and  levels 
of  from  two  to  three  times  this  value  would  not  be  unlikely. 
The  orientation  of  the  prestress  was  pure  shear  in  the 


An  approximate  relation  between  the  F  factor,  used  to 
characterize  the  anomalous,  or  tectonic,  component  of  the 
radiation  relative  to  the  direct  explosion  component  is 


(Toksoz,  et  al.,  1965): 
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(o)  DIRECT  P  WAVE  FROM  EXPLOSION  (Pg) 
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(f)  COMPOSITE  FIRST  MANTLE  ARRIVAL  SERIES 

Figure  1.  Theoretical  congressional  (P)  wave  phases  generated 
by  the  Handley  explosion  and  associated  tectonic  release.  The 
Handley  event  (1100  kt)  modeled  by  the  explosion  source  model 
T-l.  Prestress  (0^3)  for  tectonic  release  taken  to  be  65  bars, 
shatter  zone  radius  R0  =  750  meters.  Earth  structure  CIT  109- 
Low  Q  model.  Distance  4066  km,  azimuth  30°.  Vertical  component 
LRSM  short  period  seismometer,  source  depth  1.2  km. 


horizontal  plane,  so  that  relaxation  of  stress  around  the 
shatter  zone  is  equivalent  to  a  strike-slip  double  couple 
point  source.  For  a  homogeneous  prestress,  the  stress 
relaxation  around  a  spherical  shatter  zone  is  such  as  to 
always  produce  pure  quadrupole  radiation  -  that  is  a  simple 
double  couple  equivalent.  The  figure  shows  the  time  domain 
pulse  contributions  to  the  overall  P  wave  train  and  the 
important  point  is  that  the  direct  P  wave  from  tectonic 
release  is  much  smaller  than  the  explosion  generated  P  wave. 
On  the  other  hand  the  tectonic  S  wave  is  comparable  to  the 
explosion  P  wave,  but  only  contributes  upon  reflection  at 
the  free  surface  so  that  it  influences  the  wave  train  at 
later  times.  The  ratio  of  P  to  S  wave  production  by  tectonic 
release  of  this  type  scales  is: 


AP  ,Vs  3 

A  a  (V_) 
s  p 


where  A  and  A  are  the  amplitudes  of  tectonic  P  and  S  waves, 
p  s 

and  V  and  V  are  the  P  and  S  wave  velocities  in  the  source 
P  s 

region.  In  general  then,  we  expect  the  tectonic  P  wave  to 
be  of  the  order  3/1  down  in  amplitude  from  the  S  wave.  If 
we  consider  explosions  with  F  factors  near  unity,  so  that 
the  tectonic  energy  released  is  nearly  the  same  as  the 
explosive  energy  converted  to  seismic  radiation,  such  as  in 
the  example  in  Figure  1,  then  the  tectonic  P  wave  will  be  five 
times  smaller  than  the  direct  explosive  P  wave  (i.e.,  the 
energy  of  both  sources  is  about  the  same  but  80%  of  the 
tectonic  source  energy  is  contained  in  the  S  waves  produced) . 
The  example  in  Figure  1  illustrates  this  relationship.  Thus, 
as  noted  earlier,  if  the  body  wave  magnitude  measurement  is 
confined  to  the  first  cycle  of  P  wave  motion,  then  the 
effects  of  tectonic  release  will  be  minimized.  A  similar 


conclusion  was  reached  by  Bache  (1976a)  from  a  series  of 
modeling  experiments  in  which  P  waves  from  a  number  of 
explosions  were  studied,  with  tectonic  release  effects 
included  to  achieve  detailed  fits  to  the  observations. 


Displacement  Amplitude  Spectral  Density  (mlcrons/cps) 


Frequency(Hz)  Frequency(Hz) 


(a)  DIRECT  P  WAVE  FROM  (b)  REFLECTED  P  WAVE  FROM 
THE  EXPLOSION  THE  EXPLOSION 


Frequency(Hz)  Frequency  (Hz) 


(c)  SUM  OF  DIRECT  AND  ( d)  SUM  OF  DIRECT,  R  EFLECTED 
REFLECTED  P  WAVES  ANDCONVERTED  PHASES  DUE' 

DUE  TO  THE  EXPLOSION  TO  ALL  SOURCES  (explosion 

ond  tectonic  releose) 

Figure  2:  Theoretical  spectra  of  explosion  generated  compress ion al 
arrivals  from  the  Handley  event  at  a  distance  of  4066  km  and  azimuth 
of  30e .  CIT  109  Low  Q  model,  explosion  model  T-l,  source  depth  1.2  km. 
Vertical  and  radial  refer  to  the  two  components  of  the  vector  field 
at  the  earth's  surface  with  radial  denoting  the  horizontal  component 
in  the  radial  direction  from  the  so nr« 
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Figure  2  shows  the  spectra  of  the  waves  represented 
in  Figure  1,  with  the  last  inset  (d)  showing  the  spectral 
perturbation  resulting  from  tectonic  effects.  It  should  be 
emphasized  that  the  composite  spectrum  shown  in  (d)  is  for 
the  entire  P  wave  train.  If  spectral  methods  are  used,  say 
narrow  band  filtering  methods,  to  obtain  a  spectral 
magnitude  for  the  first  arrival  pulse  (i.e.,  the  first  cycle 
of  the  P  wave  train) ,  then  the  result  obtained  would  look 
like  inset  (c) ,  rather  than  (d) .  Figure  3  shows  the 
character  of  the  entire  predicted  P  wave  train,  including 
multiple  mantle  arrivals  to  be  expected  in  the  distance 
range  near  4000  km,  for  the  explosion  above  compared  with 
explosion  plus  tectonic  release.  Finally  Figure  4  shows  a 
comparison  with  an  observation,  from  the  underground  test 
Bilby.  The  predicted  seismogram  for  the  P  wave  train  is 
remarkably  similar  to  that  observed,  and  this  kind  of 
agreement  is  not  unusual.  It  is  clear  however,  that  the 
effect  of  tectonic  release  on  the  P  waves  is  not  large  and 
that  an  explosion  by  itself  could  fit  the  observations 
adequately,  especially  when  uncertainties  in  structure  and 
the  predicted  explosive  source  function  itself  are  taken  into 
account.  This  is  also  illustrated  in  Figure  3  by  the  small 
differences  between  the  theoretical  seismograms  with  and 
without  tectonic  release. 

The  predictions  of  seismic  radiation  from  tectonic 
release  have  usually  assumed  a  uniform  prestress  condition 
in  the  medium  prior  to  the  creation  of  the  explosive  shatter 
zone,  with  the  exception  of  the  treatment  employing  a  stress 
relaxation  cut-off  at  some  radius  (R  )  in  order  to  approximate 
the  effects  of  a  stress  concentration  in  the  medium  (e.g., 
Archambeau,  1970,  1972) .  When  a  uniform  prestress  (extending 
to  infinity)  is  used  and  a  spherical  shatter  zone  is  created, 
then  the  radiated  field  is  pure  quadrupole  and  its  far  field 
spectrum  is  flat  from  zero  frequency  to  a  corner  frequency, 


Figure  3:  Comparison  of  theoretical  signals  at  4066  km  from  Handley  (1100  kT)  explosion 
model  with  and  without  tectonic  release.  Explosion  model  was  T-l,  tectonic  release  based 
on  65  bar  pure  shear  prestress  (CJ13)  3nd  shock  generated  shatter  zone  of  750  meters  in 
radius,  source  depth  1.2  km. 
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Figure  4:  Comparison  of  observed  and  theoretical  seismograms 

for  “the  Bilby 'event  at  HNME,  4066  from  ground  zero 
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of  tectonic  release  (e.q.,  Archambeau  and  Sammis,  1970; 
Archambeau,  1972;  Aki  and  Tsai,  1972;  Toksoz  and  Kehrer, 

1972) .  The  only  question  seems  to  be  precisely  what  mechanism 
of  tectonic  release  is  responsible,  as  noted  earlier. 

One  line  of  evidence  concerning  the  process  involved 
in  tectonic  release  is  furnished  by  direct  comparison  of 
seismic  radiation  from  an  explosion  and  an  earthquake  occurring 
in  the  same  region.  One  such  pair,  the  Fallon  earthquake 
(m^  =  4.4,  depth  15  km)  and  the  Shoal  underground  nuclear 
explosion  test  (m^  -  4.9;  test  medium,  granite?  F  =  .58) 
was  studied  in  some  detail  by  Lambert  et  al. ,  1972.  These 
events  were  separated  by  only  a  few  tens  of  kilometers  so 
that  the  average  tectonic  environment,  as  measured  by  the 
long  period  surface  waves,  should  be  comparable.  Figure  10 
shows  one  of  the  striking  differences  between  true  earthquake 
long  period  radiation  and  the  anomalous  long  period  radiation 
from  explosions.  In  particular,  the  ratio  of  Love  to 
Favleigh  wave  (L/R)  spectral  amplitudes  strongly  increases 
with  increasing  period  for  the  earthquake,  while  the  same 
ratio  is  essentially  constant  for  explosions.  The  L/R 
spectral  ration  for  the  underground  test  Bilby  (m^  =  5.8, 
test  medium,  tuff,  F  =  .5)  is  also  shown,  for  comparison 
with  another  (much  larger)  explosion  event,  and  it  also 
shows  a  very  different  period  variation  for  L/R.  The  con¬ 
clusion  to  be  drawn  is  that  the  source  of  the  anomalous  ex¬ 
plosion  radiation,  which  is  totally  responsible  for  the  Love 
wave  excitation,  must  be  quite  different  from  a  small  earth¬ 
quake,  like  the  Fallon  event.  In  particular  it  must  be 
less  efficient  as  a  long  wave  length  radiator  than  is  a 
small  earthquake  and  thus  a  smaller  source  dimension  is  implied. 
The  tectonic  stress  relaxation  induced  by  shatter  zone  creation 
is  a  source  of  relatively  small  characteristic  dimension  that 
produces  quadrupole  long  period  radiation  with  considerable 
SH  wave  production.  Hence,  it  is  a  most  likely  model  for  the 


Elastic  wave  radiation 


Fig.  /ft.  The  variation  of  L/R  with  period  for  the  BILBY  and  SHOAL  explosions 
and  the  Fallon  earthquake,  in  the  distance  range  near  800  km.  L/R  for  the 
earthquake  is  larger  than  2-0  and  increases  with  increasing  period.  This  difference 
in  behaviour  is  a  consequence  of  the  much  larger  source  dimension  associated 
with  the  earthquake. 


Fro.  ft  Radiation  patterns  of  Love  and  Rayleigh  waves  from  BILBY  at  periods 
of  15  and  20  s.  The  insets  show  the  theoretical  pattern  shapes  obtained  as  a  fit  to 
the  observations,  using  superposed  monopole  and  quadrupole  point  sources 
with  fixed  relative  excitation  of  quadrupole  to  monopole  of  F  =  0-5,  with 
quadrupole  1  strike '  »  342“.  The  point  source  equivalents  are  the  same  for 

both  Love  and  Rayleigh  waves.  The  factor  e  =  0-7  is  the  particle  orbit  ellipticity 
factor  for  the  Rayleigh  waves,  which  depends  on  the  structure  used  in  the 
theoretical  calculations. 
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anomalous  explosion  source  component.  Figure  11  shows  the 
form  of  the  observed  radiation  patterns  for  Love  and  Rayleigh 
waves  at  periods  T  =  15  sec  and  20  sec.  The  insets  are  the 
theoretically  predicted  surface  wave  patterns  based  on  tectonic 
release  due  to  the  explosion  produced  spherical  shatter 
zone,  plus  the  pure  explosion  monopole  field  (Lambert  et  al., 
1972) .  The  prestress  orientation  used  is  such  as  to  be 
consistent  with  the  predominently  strike  slip  mechanism 
associated  with  the  Fallon  earthquake,  while  the  magnitude  of 
the  initial  stress  required  to  fit  the  amplitude  of  both 
the  Rayleigh  and  Love  waves  was  65  bars  (Archambeau  and 
Sammis,  1970).  The  (mean)  prestress  level  seems  entirely 
reasonable.  Thus  the  "shatter  zone  model"  seems  consistent 
with  the  observations  of  surface  wave  radiation,  in  that  the 
values  of  prestress  magnitude  and  orientation  required  to 
fit  both  Love  and  Rayleigh  wave  observations  are  compatable. 
with  the  tectonics  of  the  region. 

Figures  12  and  13  show  how  well  this  model  simultaneously 
fits  the  long  period  surface  wave  radiation  for  the  Shoal 
and  Bilby  explosions.  Here  observed  L/R  ratios  at  several 
stations  (where  the  maximum  observed  surface  wave  amplitudes 
are  near  15  seconds)  are  compared  to  the  predicted  ratio  at 
different  azimuths.  The  agreement  is  very  good,  considering 
the  probable  lateral  refraction  effections  to  be  expected. 
Similar  analysis  for  numerous  other  explosions,  for  example 
by  Toksoz  and  Kehrer,  1972,  shows  similar  results. 

Hence  it  would  appear  that  the  production  of  seismic 
radiation  by  spherical  shatter  zone  induced  (tectonic)  stress 
relaxation  can  explain  the  anomalous  long  period  surface 
wave  observations,  and  the  P  and  S  wave  train  complications 
as  well,  in  the  great  majority  of  observed  cases.  Particularly 
for  the  long  period  surface  waves,  which  are  relatively 
unaffected  by  small  dimensional  high  stress  concentrations, 
it  should  be  possible  to  estimate  the  magnitude  of  the 
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Fig.  /£.  The  theoretical  L/R  ratio  for  SHOAL,  at  IJ-s  period,  as  a  function  of 
azimuth  and  using  the  source  parameters  ni^n^c-d-  Observations  at  various 
azimuths  where  signal-to-noise  ratios  where  high  are  indicated  and  identified  by 
abbreviated  station  symbols.  This  also  shows  the  nature  of  the  L/R  azimuth 
variation  (which  is  nearly  identical  to  that  for  BILBY)  and  that  a  single  equivalent 
source  fits  both  Love  and  Rayleigh  wave  radiation  for  this  event. 


Fig./S,  The  theoretical  L/R  for  BILBY,  at  T  *»  15  s,  as  a  function  of  azimuth 
using  the  source  parameters  i/jjitAfeJ- Fhe  observed  ratio  at  a  number  of  stations 
with  good  signal-to-noise  ratios  Tot  both  Love  and  Rayleigh  waves  are  shown  at 
their  appropriate  azimuths.  This  also  shows  the  nature  of  the  L/R  azimuth 
variation,  and  it  is  clear  that  a  single  equivalent  source  fils  both  Love  and  Rayleigh 
wave  radiation  for  this  event. 
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anomalous  effects  (perturbations)  to  be  expected.  On  the 

other  hand  the  body  wave  perturbations  are  much  more  difficult 

to  predict,  because  of  their  dependence  on  details  of  the 

spatial  dependence  of  the  prestress,  but  the  large  effects 

can  be  avoided  if  only  the  first  cycle  of  the  P  wave  train 

* 

is  used  in  the  event  analysis. 

Figure  14  (Archambeau  et  al.,  1974)  shows,  among 
other  things,  the  effects  of  tectonic  stress  on  mb  and  Mg 
for  explosions  in  granite,  covering  the  yield  range  from  1  KT 
to  1000  KT,  where  shatter  zone  stress  relaxation  is  assumed. 

The  propagation  path  is  appropriate  to  NTS,  so  that  the  upper 
mantle  includes  a  well  developed  low  velocity  and  low  Q 
zone.  For  this  reason  the  theoretical  mb  values  are  some¬ 
what  lower  than  would  be  expected  in  a  non- tectonic  region, 
but  are,  on  the  average,  what  appear  to  be  appropriate  for 

NTS  and  the  Basin  and  Range  Geologic  Province.  The  circles 

*  * 

represent  predictions  of  mb  and  Mg  values  from  the  explosions 
alone,  without  tectonic  release.  The  measurements  at  these 
points  are  made  from  the  synthetic  seismograms,  using  the 
first  cycle  of  the  P  wave  train  (denoted  m^^  )  and  the 
Rayleigh  wave  vertical  component  at  20  seconds.  The  "error 
bars"  (upward  vertical  lines  and  horizontal  lines  attached  to 
each  circle)  indicate  the  possible  increase  in  Mg  if  the  Airy 
phase  is  measured  and,  for  m^,  the  increase  to  be  expected 


★ 

The  situation  for  the  prediction  of  anomalous  effects  in  this 
context  is  similar  to  the  situation  described  by  Bismark 
concerning  Prussia  and  Austria;  to  paraphrase:  "The  situation 
(for  long  period  surface  waves)  is  serious  but  not  hopeless, 
the  situation  (for  short  period  body  waves)  is  hopeless  but 
not  serious." 

•k  *  .  , 

The  superscript  "R"  is  used  on  Ms  in  the  figure  to  identify 
the  Ms  measurement  as  being  associated  with  the  Rayleigh  wave. 
This  was  employed  because  a  superscript  "L"  was  used  for  the 
similar  magnitude  measurement  from  Love  waves.  This  notation 
is  unnecessary  here  as  Ms  will  always  be  the  standard  magnitude 
from  the  Rayleigh  wave. 


mb 


Fig.  14.  Rayleigh  surface  wave  magnitude  as  a  function  of  body 
wave  magnitude  for  explosions  in  granite  (model  G-l) . 

The  circles  denote  theoretical  values  of  m^l)  and  M 
for  zero  prestress,  the  squares  for  165  bar  prestrels. 

The  horizontal  lines  indicate  the  range  of  possible  mb 
values,  the  vertical  lines  indicate  the  range  of  possible 
Ms  values.  The  large  variation  in  Ms  arises  from  the 
possibility  of  measuring  an  Airy  phase.  The  variation 
in  mb  arises  from  the  possibility  of  measuring  the  2nd 
or  3  rd  cycle  in  the  P  wave  train  rather  than  the  first. 
The  upper  and  lower  bound  lines  for  earthquake  data  are 
shown  for  comparison  along  with  the  mean  NTS  explosion 
magnitude  line.  The  X  symbols  are  individual  NTS 
explosions  defining  the  extreme  upper  range  of  observed 
M  values  for  explosions  from  Nevada.  The  yields  of 
the  theoretical  explosions  are  indicated  along  the  m. 
scale.  b 
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if  the  larger  of  the  second  or  third  cycle  in  the  P  wave 

train  is  used  to  measure  The  squares  correspond  to  the 

and  Mg  values  that  would  be  obtained  with  165  bars  prestress. 

The  Mg  value  here  is  not,  however,  an  average  over  the  whole 

radiation  pattern,  but  the  largest  that  could  be  measured. 

It  is  therefore  an  upper  bound  on  the  (positive)  perturbation 

in  Ms  that  can  occur  for  this  kind  of  tectonic  stress  release. 

Naturally  the  "error  bars"  are  dragged  along  with  the  point 

in  the  m^  -  Mg  plane  and  the  upper  horizontal  bar  is  associated 

with  the  limit  Mg  measured  from  the  Airy  phase.  The  horizontal 

extent  of  the  bar  indicates  the  perturbation  possible  in 

due  to  "cycle  selection".  The  prestress  orientation  was 

taken  to  be  such  that  the  radiation  quadrupole  corresponded 

to  a  "strike-slip"  equivalent  double  couple.  For  other 

choices  of  prestress  orientation, the  maximum,  and  average 

as  well,  could  be  much  reduced  instead  of  being  increased. 

(e.g.,  A  thrust  equivalent  could  probably  cancel  or  even 

reverse  the  directly  generated  explosion  Rayleigh  waves.) 

In  any  case  the  figure  illustrates  the  size  of  the  effects 

for  a  particular  case,  and  they  are  significant  for  M  . 

They  are  not  very  large  for  m^  and  indeed  the  variations  in 

don't  change  much  with  prestress  in  this  special  case. 

However , the  calculation  was  done  with  essentially  uniform 

prestress  and  when  stress  concentrations  are  present,  as  they 

almost  certainly  would  be  in  actuality,  the  m.  variations 

could  be  much  larger.  Note  also  that  the  m^  values  are 

essentially  identical  with  and  without  prestress,  as  shown 

in  the  figure.  (i.e.,  The  circle  and  square  points  are 

only  displaced  vertically,  showing  only  a  perturbation  in 

M  at  20  seconds.)  While  there  actually  were  slight  changes 
s  ( 1 ) 

in  the  m^  values  for  the  two  cases,  they  were  too  small 
to  be  shown  on  this  plot.  However,  with  stress  concentrations 
present  they  could  be  considerably  larger. 
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Also  shown  on  the  figure  are  lines  indicating  earthquake 
populations  for  different  regions,  and  points  (x)  denoting 
some  of  the  more  anomalous  NTS  explosions  (generally  those  hav¬ 
ing  high  Ms  values  relative  to  their  ro^  values) .  Note  that  some 
of  these  events  fall  close  to  the  extreme  values  for  Ms  pertur¬ 
bations  due  to  this  kind  of  tectonic  release.  Comparing  them 
with  the  earthquake  population  limit  lines  shows  that  they  could 
be  confused  with  earthquakes  on  the  basis  of  the  m^-Mg  discrim¬ 
ination  criteria  alone.  Finally,  for  general  reference,  the 
mean  line  for  explosions  at  NTS  is  also  shown. 


Summary:  State  of  Knowledge  and  Research  Needs 

The  previous  discussion  is,  in  effect,  a  summary  of  what 
is  considered  to  be  the  current  state  of  knowledge.  In  addition, 
it  is  to  be  hoped  that  the  uncertainties  and  ambiguities  are 
reasonably  well  covered.  By  way  of  a  summary,  then,  if  seems 
most  useful  to  briefly  state  the  essential  conclusions  to  be 
drawn,  albeit  with  some  being  rather  tenative,  and  to  then  list 
areas  of  research  that  could  provide  the  necessary  details  for 
very  firm  conclusions. 

The  principal  conclusions  and  results  are: 

(1)  Stress  relaxation  effects  can  be  expected  in 
any  material  capable  of  sustaining  long  term  non¬ 
hydrostatic  stress.  The  largest  effects  will  oc¬ 
cur  in  the  (known)  regions  of  high  tectonic 
activity.  That  is,  large  effects  would  be  expected 
at  plate  margins  in  and  near  intrusive  zones  and 
generally  where  loading  of  the  crust  is  evident, 
such  as  near  large  river  deltas.  Because  of  the 
likelihood  of  high  stresses  in  seismicallv  quiet 
zones  along  plate  boundaries,  it  is  probable  that 
seismic  gaps  are  areas  of  special  importance. 

(2)  It  appears  that  most  of  the  well  documented 
anomalous  effects  in  the  seismic  radiation  field 
from  explosions  (large  SH  wave  production,  observed 
especially  at  long  periods  from  10-25  seconds; 
strong  perturbations  in  the  long  period  Rayleigh 
type  surface  waves,  azimuthally  dependent  increase 
in  complexity  of  the  short  period  P  waves)  are  due 
to  tectonic  stress  relaxation  in  the  vicinity  of 


the  roughly  spherical  shatter  zone  created  by  the 
explosive  shock  wave.  However,  there  may  be,  at 
least  in  some  cases,  strong  asymmetries  in  the 
fracture  zone  around  an  explosion;  especially  when 
the  medium  is  highly  stressed  and/or  existing  stress 
loaded  faults  are  nearby.  In  this  case  triggering 
of  an  earthquake  is  said  to  have  occurred.  However, 
the  likdihood  of  this  occurrence  is  judged  to  be 
small  and  to  have  rarely  occurred  in  the  past  on 
a  large  scale.  Nevertheless,  it  is  possible  that 
future  tests  could  be  so  arranged  so  as  to  maximize 
the  likelihood  of  such  an  occurrence.  When  an  earth 
quake  at  or  very  near  the  test  point  is  induced, 
discrimination  and  especially  yield  estimation, 
would  be  considerably  more  difficult. 

(3)  Long  period  surface  wave  radiation  is  strongly 
perturbed  by  tectonic  release  effects  within  the 
whole  measureable  low  frequency  band  (i.e.,  from 
approximately  5  sec.  to  60  sec.  in  period).  The 
perturbations  in  the  observed  Rayleigh  wave  forms 
can  be  such  as  to  add  to,  or  subtract  from,  the  ex¬ 
plosive  generated  Rayleigh  wave  depending  on  the 
orientation  and  magnitude  of  the  prestress  in  the 
vicinity  of  the  explosion.  The  magnitude  of  this 
effect  can  be  very  (unacceptably)  large.  In  those 
cases  where  Love  waves  are  significant,  so  that 
tectonic  release  is  involved,  then 

yield  estimation  using  Ms  can  only  be  made  after 
correction  of  the  Rayleigh  wave  measurement  using 
the  observed  Love  wave  to  deduce  the  size  and  con¬ 
figuration  of  the  tectonic  source.  Such  a  correc¬ 
tion  would  be  much  more  reliable, when  spherical 
shatter  zone  induced  tectonic  release  is  involved, 
than  it  would  be  if  actual  earthquake  triggering 
is  involved. 

(4)  Short  period  perturbations  in  the  wave  train 
can  be  expected  to  be  very  complex  due  to  dependence 
on  stress  concentration  effects  and  local  complex 
structure.  However,  the  perturbations  should  be 
small  to  moderate  for  the  first  cycle  of  the  P 

wave  motion,  while  being  significantly  larger  for 
the  later  part  of  the  P  wave  train.  Body  wave  mag¬ 
nitude  measured  from  the  first  P  wave  cycle  should, 
therefore,  be  minimally  perturbed  by  stress  relaxa¬ 
tion  effects  and  the  complexities  of  local  structure 
Further,  corrections  to  the  first  cycle  of  the  P 
wave  for  tectonic  affects  could  conceivably  be  made 
for  purposes  of  yield  estimation. 
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Appendix  2 


Regional  Attenuation  Effects  on  P  Waves  and  Effects 
of  Attenuation  on  Surface  Waves 
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Regional  Attenuation  Effects  on  P  Waves  and  Effects  of  Attenuation  on  Surface 
Waves .  Charles  Archambeau 

Introduction. 

The  effects  of  attenuation  for  both  surface  and  body  waves  can  be  described, 
most  appropriately,  in  terms  of  intrinsic  dissipation  functions  Q0(r,  f) 
and  Qc(r,  f)  for  the  earth.  Here  Q  and  Q0  are  Q  functions  for 

P  Ct  p 

compressional  and  shear  wave  losses  respectively,  and  both  are  functions  of 
radius  and  frequency+.  They  are  also  functions  of  the  other  spatial  coordinates, 
but  it  is  easier  and  actually  most  appropriate  to  define  different  Q  models 
for  different  geologic  provinces.  It  is  quite  clear  from  observational 
results  for  body  waves  (e.g.,  Archambeau,  Flinn  and  Lambert,  1969)  that  much 
of  the  dissipation  of  P  waves  takes  place  in  the  low  velocity  zone. 

Therefore  high  attenuation  is  correlated  with  high  heat  flow  and  large  P 
delays,  these  in  turn  characterizing  geophysical-geological  provinces. 

In  particular,  shield  areas  with  low  heat  flow  and  negative  P-delays  show 
low  attenuation,  while  active  tectonic  provinces  with  high  heat  flow  and 
positive  P-delays  show  high  attenuation.  All  of  these  effects  are  clearly 
related  to  the  depth  span  and  intensity  of  the  low  velocity  zone.  These 
correlations  are  clearly  shown  form  the  studies  of  pP  pulses  from  earth¬ 
quakes  in  trench  zones  (e.g.,  Barazangi,  Pennington  and  Isacks,  1975), 
as  well  as  from  teleseismic  P  wave  observations  from  explosions  (e.g., 

+Qa  and  Qg  can  be  related  to  each  other  under  the  assumption  that 
dissipation  in  pure  compression  is  very  small  relative  to  losses  in 
shear.  Then,  for  typical  mantle  elastic  velocities,  Q  =  9/4  Q„  . 

For  details  see  Anderson  et  al.,  1965.  P 


Der  et  al.,  1975;  Der  and  McElfresh,  1976).  Thus,  for  body  wave  magnitudes, 
one  expects  variations  in  which  are  directly  related  to  the  geologic 
provinces  of  the  source  and  receiver.  In  particular,  sources  in  tectonic 
provinces  will  show  reduced  m^  values  at  1  Hz  relative  to  the  same  sources 
in  shield  regions.  Since  tectonic  provinces  show,  in  general,  highly 
variable  low  velocity  zone  thicknesses  and  correspondingly  variable 
heat  flow  values  and  P  delays,  one  can  also  expect  variability  in  the  m^ 
reduction,  from  quite  large  reductions  to  rather  small  reductions,  depending 
on  precisely  where  the  source  is  located. 

Similar  statements  can  be  made  about  surface  wave  attenuation  relative 
to  Geological-Geophysical  provinces.  That  is,  the  strongest  attenuation 
occurs  within  the  low  velocity  zone  and  strong  surface  wave  attentuation 
is  correlated  with  regions  of  high  heat  flow,  large  P  (and  S)  wave  delays 
and  tectonic  activity.  For  shorter  period  surface  waves  not  penetrating 
the  low  velocity  zone  however  (i.e.,  for  periods  less  than  30  seconds), 
the  attenuation  is  not  as  large  as  for  longer  period  surface  waves  (i.e., 
the  observed  Q  is  about  300  compared  to  observed  Q  values  of  around 
100  for  surface  waves  in  the  period  range  30-200  seconds)  and  there  is 
less  regional  dependence  in  attenuation  (e.g.,  Solomon,  1972).  However, 
Mitchell,  1975,  has  shown  that  for  rather  short  period  surface  waves,  near 
5  seconds,  the  attenuation  is  quite  strongly  regionally  dependent.  Never¬ 
theless,  he  finds  that  for  the  longer  periods  up  to  30  seconds,  there  is 
little  regional  dependence.  This  very  short  period  regional  dependence 
in  attenuation  is  probably  more  related  to  scattering  than  to  anelastic 
effects,  in  that  tectonically  active  provinces  usually  show  larger  near  surface 
lateral  variability  in  velocity  structure  than  do  the  more  stable  provinces. 
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Thus  M  values  based  on  20  second  Rayleigh  waves  do  not  show  strong  regional 
s 

"Q-bias". 

The  frequency  dependence  of  the  anelastic  dissipation  has  only  recently 

been  considered  in  any  great  detail.  Originally  Archambeau  et  al.,  1969, 

showed  that  P^  phases  in  the  Western  U.S.  were  attenuated  such  that  the 

high  frequencies  required  higher values  than  the  lower  frequencies  - 

that  is,  the  Qq  appeared  to  increase  with  frequency  in  the  frequency  range 

from  .5  to  3  Hz.  These  direct  observations  were  also  in  agreement 

with  the  observation  that  Q  (and  Q0  )  models  obtained  from  low  frequency 

ot  p 

surface  waves  had  lower  Qq  values,  essentially  everywhere  in  the  mantle, 
when  compared  to  the  model  obtained  from  high  frequency  (1  to  3  Hz)  body 

wave  observations.  The  upper  mantle  Qq  models  that  have  been  obtained 
from  low  frequency  surfa-.  wave  and  free  oscillation  data  and  from  high 
frequency  body  wave  data  are  shown  in  Figure  1.  The  model  SL8  is  from 
the  analysis  of  free  oscillation  data  by  Anderson  and  Hart  ,  1978;  the 
model  MM8  is  from  surface  wave  data  inversion  by  Anderson  et  al.,  1965; 
and  the  model  AFL  is  from  body  wave  data  inversion  by  Archambeau  et  al., 

1969.  Each  model  applies  only  to  the  frequency  range  covered  by  the  data  used 
to  obtain  it.  The  trend  of  these  results  is  toward  high  Q  values  with 
increasing  frequency  of  the  data  used  in  this  inversion. 

Solomon,  1972b,  proposed  a  frequency  dependent  intrinsic  Q  for  the 
mantle  involving  activated  processes  that  satisfied  the  observed  Long  period 
surface  wave  dispersion  quite  well.  Liu  and  Archambeau,  1975  and  1976, 
showed  that  this  model  fit  the  total  set  of  surface  wave  and  free  oscillation 
data  quite  well  and  that  it  predicted  relatively  large  shifts  in  the 
dispersion  (group  velocity  versus  frequency)  and  free  oscillation  periods. 
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showing  that  the  effective  velocity  structure  sensed  by  low  frequency  waves 
is  different,  than  that  for  higher  frequency  waves.  Liu  et  al.»  1976, 
expanded  upon  these  results  and  proposed  an  absorption  band  intrinsic 
Q  model  that  consisted  of  a  distribution  of  activated  processes,  each  with 
a  different  characteristic  relaxation  time  corresponding  to  a  superposition 
of  many  absorption  processes  acting  to  dissipate  energy.  This  model  was 
also  shown  to  be  compatible  with  observations  of  surface  wave  and  free 
oscillation  dissipation.  The  absorption  and  Q  model  amounts  to  an  ex¬ 
tension  of  Solomon's  model,  wherein  many  activated  processes  are  allowed 
racher  than  one  or  two,  and  is  more  realistic  in  terms  of  the  known  micro¬ 
physics  of  crust-mantle  materials. 

Currently  this  kind  of  intrinsic  Q  model  is  being  used  to  constrain 
the  frequency  dependence  of  the  intrinsic  Q  in  the  earth,  in  order  to 
invert  for  both  the  depth  dependence  and  the  shape  of  the  absorption  band 
at  each  depth  (and  hence  the  intrinsic  frequency  dependence  of  Qa  and 
Qd  ).  Figure  2,  from  Lundquist,  1980,  shows  the  form  of  the  absorption 

D 

band  models  being  used.  Such  an  absorption  band  applies  at  each  depth  in 
the  earth  and  varies  with  temperature  pressure,  material  chemistry  and 

phase  state.  The  parameters  and  are  low  and  high  frequency 

"relaxation  times"  corresponding  to  the  half  amplitude  points  on  the 
"Q-filter"  in  the  frequency  domain.  These  parameters  are  treated  as  unknowns 
and  are  obtained,  as  functions  of  depth,  by  inverting  the  observed  attenuation 
data. 

The  Q  models  shown  in  Figute  1  have  very  poor  resolution  of  Q  variations 


in  the  crust,  mainly  because  little  or  no  very  short  period  surface  wave  data 
was  used  for  the  inversion  with  the  surface  wave  and  free  oscillation  data 
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and,  in  the  case  of  the  body  wave  derived  model,  the  sampling  of  the  crust 

using  teleseismic  P  waves  was  minimal.  Mitchell,  1980,  has  however  studied 

relatively  short  period  surface  wave  propagation  in  the  Eastern  U.S.  (with 

periods  from  1  to  40  seconds)  and  obtained  crustal  Q0  models  in  some 

P 

detail.  He  has  also  shown  that  the  intrinsic  Q„  can  be  best  described 

P 

by  a  frequency  dependent  Q-f unction  of  the  form: 

Qg(w,  r)  *=  Qg(r)0JC 

with  X,  between  .3  and  .5  for  the  period  range  1  to  40  seconds.  Here  again 
it  appears  that  the  intrinsic  Q  increases  with  increasing  frequency,  however 
such  a  conclusion  based  on  the  fits  given  by  Mitchell  may  be  premature. 

In  any  case  his  models  show  a  Qg  average  of  about  250  in  the  upper  crust 
(0-15  km)  and  near  1000  for  the  lower  crust  (15-40  km).  These  values  are 
significantly  higher  than  the  Qg  values  in  the  low  velocity  zone  of  the 
upper  mantle,  where  Qg  -  50-100  is  appropriate. 

It  is  clear  that  the  effects  of  attenuation  on  surface  wave  magnitudes, 
measured  at  20  seconds,  are  not  as  extreme  as  are  attenuation  effects  on 
body  wave  magnitudes.  First,  there  is  little  observed  regional  variation 
in  attenuation  in  this  period  range.  Second,  the  attenuation  is  not  very 
large,  that  is  the  Q  of  the  crust,  while  of  course  variable  in  both 
frequency  and  with  depth,  is  quite  high.  Thus,  corrections  in  Mc  for 
attenuation  could  be  made  and  they  would  not  be  very  large.  It  is  of 
course  important  that  Mg  be  measured  at  20  seconds. 

Frequency  Dependent  Q  Models  for  Teleseismic  P  waves  and 
Mantle  Surface  Waves 

The  best  (i.e.  only)  first  order  frequency  dependent  Q  model  for 
the  upper  mantle  has  been  obtained  by  Lundquist,  1980,  The  model  uses 


an  absorption  band  intrinsic  Q  of  the  type  shown  in  Figure  2.  The  model 

is  obtained  by  first  taking  the  previously  determined  low  frequency  Q  models 

(the  MM8  and  SL8  models  in  Figure  1)  and  the  high  frequency  model  AFL  in 

Figure  1  as  appropriate  Q  variations  in  the  mantle  in  the  frequency  ranges 

for  which  they  are  defined.  That  is,  the  frequency  dependent 

model  is  constrained  to  give,  to  first  order  at  least,  the  SL8  model  at 

very  low  frequencies  and  the  AFL  model  at  high  frequencies,  near  3  Hz. 

The  observed  Q  models  in  Figure  1  turn  out  to  imply  that  there  is  one 
absorption  band  model  for  the  mantle  beneath  the  low  velocity  zone,  having 
regular  properties  varying  with  depth  in  a  manner  consistent  with  the 
temperature-pressure  variations  in  the  earth  in  this  depth  range,  and  a 
separate,  very  different,  kind  of  absorption  band  which  appears  to  be 
confined  to  the  low  velocity  zone.  The  absorption  band  for  the  low 
velocity  zone  appears  to  be  narrow  (i.e.,  x^  and  X2  relatively  close  in 
value)  while  the  lower  mantle  absorption  band  appears  to  be  very  broad 
(i.e.,  Tj,  very  large  and  x2  near  .1  sec).  The  second  absorption 
band  associated  with  the  low  velocity  zone  may  be  a  consequence  of  a 
partial  melt  state  within  the  zone.  In  any  case  it  is  confined  to  this 
zone  and  therefore  varies  with  the  extent  and  intensity  of  the  low  velocity 
zone. 

Using  such  a  rough  double  absorption  band  Q  model  as  a  starting  point, 
the  frequency  dependent  Q  model  can  be  refined  by  adjusting  the  various 
absorption  band  parameters  (in  particular  the  "high  and  low"  frequency 
relaxation  times  x2  and  x^  plus  the  maximum  Q  ^  level  of  the  absorption 
band  at  each  depth)  to  fit  frequency  and  time  domain  observations.  In  particular, 
Lundquist  adjusted  the  starting  double  absorption  band  model  to  be  such  that 
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when  attenuation  corrections  are  applied  to  observed  earthquake  and  explosion 
P  wave  spectra,  then  corrected  source  spectra  had  high  frequency  asymptotic 
behavior  of  the  form  1/oj  or  l/o a  ,  as  is  expected  from  source  theory 
considerations.  Further,  he  used  the  resulting,  somewhat  refined,  Q  model 
to  predict  time  domain  synthetic  P  wave  forms  and  further  adjusted  the 
model  to  achieve  detailed  fits  to  the  first  cycle  of  the  P  wave  train  from 
explosions.  (Only  the  first  cycle  of  the  P  wave  train  is  reasonably  well 
predicted  by  currect  explosions  models.  Further,  it  is  relatively  free 
from  uncertainties  introduced  by  near  source  structure,  tectonic  release 
and  spall  phase  production.) 

The  net  result  was  that  the  initial  double  absorption  band  model,  in¬ 
ferred  from  the  low  and  high  frequency  Q  models  of  Figure  1,  fit  the  observations 
from  NTS  explosion  events  very  well,  with  little  adjustment  necessary. 

Thus  this  model  closely  corresponds  to  the  free  oscillation  model  SL8  at  low 

frequencies  and  the  body  wave  model  AFL  at  high  frequencies  and  predicts 

the  behavior  of  a  mantle  Q  at  other  frequencies  such  that  both  spectral  and 
time  domain  observations  are  well  satisfied.  Figure  3  shows  the  properties 
of  this  double  absorption  band  model  (solid  line)  as  a  function  of  depth  and 
frequency  in  the  earth.  The  upper  inset  indicates  the  apparent  Qq  in  the 
crust,  which  is  poorly  resolved  but  is  high,  as  indicated.  The  next  inset 

shows  the  typical  form  of  the  double  absorption  band  in  the  low  velocity 

zone.  The  dotted  line  shows  the  single  absoprtion  band  that  would  exist  if 
the  low  velocity  zone  were  absent,  so  that  the  departure  of  the  solid  line 
from  the  single  absorption  band  Qq  indicates  the  effect  of  the  second 
absorption  band  associated  with  the  low  velocity  zone.  At  greater  depths  the 
variation  in  the  absorption  band  is  such  that  the  maximum  level  of  Qq  for 
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Figure  3.  Absorption  bands  in  the  Earth- vs- frequency 
and  depth.  The  dashed  line  represents  the  single 
absorption  band  model,  TF3S8,  and  the  solid  line 
represents  the  double  absorption  band  model,  TFJ12. 

At  low  frequencies,  both  models  converge  to  Q“*(SL8), 
and  model  TFJ12  converges  to  Q_1(AFL)  at  3  Hz. 
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the  band  Increases  (1/Qq  decreases)  and  the  relaxation  time  increases; 

both  uniformly  in  a  manner  controlled  by  the  temperature-pressure  increases 
with  depth. 

On  the  other  hand,  for  non-tectonic  regions, the  Qq  variation  with 
depth  and  frequency  was  found  to  be  somewhat  different  than  for  the  Basin 
and  Range  region.  In  particular,  using  explosion  data  from  Novaya  Zemlya, 
so  that  a  stable  platform  region  was  sampled,  the  Qq  variation  with 
depth  was  intermediate  between  the  single  absorption  band  variation  shown 
by  the  dotted  line  in  Figure  3  and  the  double  absorption  band  model  for  the 
Basin  and  Range.  (See  Lundquist,  1980  for  details.)  This  appears  to  be 
due  to  a  less  intense  and  thinner  low  velocity  zone  for  the  stable  platform 
region  and  a  correspondingly  more  depth  confined  and  less  intensive  second 
absorption  band  in  the  45-200  km. depth  range.  This  of  course  again  implies 
regional  variations  in  attenuation,  but  specifically  that  this  variation  is 
controlled  by  the  presence  or  absence  of  the  second  absorption  band.  Further, 
because  of  the  nature  of  this  absorption  band,  in  particular  its  frequency 
band  width,  the  frequency  dependence  of  the  absorption  can  be  quite  different 
from  region  to  region. 

k 

The  consequences  of  this  kind  of  Q  model,  in  terms  of  t  (total 

travel  time  divided  by  the  effective  Q  over  the  path  of  the  wave) ,  are 

k 

shown  in  Figure  4  for  the  double  absorption  band  model.  Clearly  "t  is 
quite  strongly  frequency  dependent. 

These  results  have  a  number  of  important  implications.  First  it  seems 

k 

evident  that  1  should  not  be  used  in  modeling  work,  but  rather  the  Qq 
or  Qg  models  should  be  used  and  modeling  should  be  done  in  the  frequency 
domain  in  order  to  properly  account  for  both  the  depth  and  frequency  dependence 
of  the  Q  and  for  the  different  apparent  "elastic"  velocities  sensed  by 
waves  of  different  frequency.  Second,  high  frequency  seismic  energy  is 
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Figure  4.  t*-vs-distance  for  selected  frequency. 
The  velocity  model  used  is  QM 2,  and  the  Q  model. is 
the  double  absorption  band  model,  TFJ12.  t*(TFJ12) 
converges  to  t*(AFL)  at  3  Hz  and  converges  to 
t*(SL8)  as  frequency  decreases. 
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propagated  with  much  less  attenuation  than  was  previously  supposed  (by  some 
at  least).  Finally,  the  efficiency  of  the  high  frequency  propagation  may 
be  highly  variable  from  province  to  province  and  will  be  correlated  with  the 
extent  and  intensity  of  the  low  velocity  zone. 

Effects  of  Attenuation  on  Surface  Waves 

The  Q  models  derived  from  body  waves  can  be  used  to  predict  the 
attenuation  of  surface  waves  and  vice-versa.  Obviously  combinations  of 
surface  and  body  wave  data  can  be  used  to  infer  Q  ^  and  Qg  as  well, and 
either  can  be  predicted  from  the  results.  Therefore  the  Qq  models 
discussed  in  the  previous  section  can  be  used  to  infer  Qg  (from  a  relation 
such  as  -  9/4  Qg  ),  and  the  resulting  model  can  be  used  to  predict,  to 
first  order  at  least,  the  expected  surface  wave  attenuation.  Lundquist's 
model  is  ,  in  this  regard,  adequate  for  the  prediction  of  the  longer  period 
surface  wave  attenuation  (T  >  40  sec)  but  is  not  well  enough  defined  in  the 
crust  to  give  very  accurate  predictions  for  shorter  period  surface  waves. 

Inasmuch  as  the  20  second  period  fundamental  model  Rayleigh  wave,  in  particular, 
is  of  major  interest  in  view  of  its  use  in  Ms  calculations,  it  is  necessary 
to  consider  high  resolution  crustal  Q  models,  such  as  are  being  obtained 
by  Mitchell  (1980).  It  seems  sufficient  here  to  only  refer  to  Mitchell's  work 
and  to  recall  the  general  comments  made  earlier  in  the  introduction.  In 
particular,  that  a  slight  frequency  dependence  is  inferred  in  Qg,with  the 
Qg  increasing  with  increasing  frequency;  that  the  mean  Qg  in  the  crust  is 
relatively  highjand  that  regional  dependence  of  the  attenuation  of  surface 
waves  in  the  range  5  <  T  <  40  sec  is  small.  All  in  all  it  does  not  appear 
that  crustal  surface  wave  attenuation  is  particularly  difficult  to  deal  with 
for  purposes  of  Mg  corrections,  to  obtain  yield  estimates  and  for  discrimination. 
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More  difficult,  and  much  larger,  are  corrections  in  Mg  for  tectonic  effects 
and  probably  for  lateral  variations  in  structure.  Uncertainties  in  these 
latter  effects  completely  overwhelm  any  correction  uncertainties  due  to 
anelasticity  effects. 

Summary:  State  of  Knowledge  and  Research  Needs. 

The  essential  conclusions  of  this  report  are: 

(1)  Frequency  dependent  Q  models  appear  to  be  required  in  both  the 
crust  and  upper  mantle.  Absorption  band  models,  with  the  Q  magnitude 
and  frequency  dependence  varying  with  depth  appear  to  be  physically 
realistic  and  to  satisfy  the  available  data. 

(2)  A  single  absorption  band  appears  to  be  appropriate  for  the 
entire  mantle  exclusive  of  the  low  velocity  zone.  Within  the  low 
velocity  zone,  when  present,  a  second  narrow  absorption  band  appears  to 
existand  accounts  for  the  increased  attenuation  and  different  frequency 
dependence  of  the  attenuation  in  tectonic  regions.  This  second  absorption 
band  is  the  likely  mechanism  for  variablility  of  body  wave  absorption 
from  region  to  region. 

(3)  High  frequency  seismic  body  waves  propagate  with  relativly 
great  efficiency  from  (and  within)  regions  not  having  a  well  developed 
low  velocity  zone.  Tectonic  zones  will  typically  absorb  much  more  of 
the  high  frequency  energy  and  this  will  generally  result  in  lower  m^ 
values.  For  this  reason  m^  should  be  measured  "spectrally"  (i.e.,  by 
narrow  band  filtering  at  1  Hz  with  the  first  cycle  of  the  P  wave  train 
selected) to  avoid  measuring  m^  at  different  effective  periods. and  Q 
corrections  should  be  made  in  order  to  account  for  differences  in  the 
regional  Q  structure. 

(4)  Surface  waves  in  the  5-40  sec  period  range  are  not  attenuated 
strongly  by  the  crustal  Q  and  there  is  no  strong  regional  dependence 

in  the  attenuation  in  this  period  range.  For  longer  periods  there  would. 
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however,  be  some  fairly  significant  regional  variations  due  to  the  variations 
in  the  low  velocity  zone.  Because  of  the  inferred  high  Q  of  the  crust, 
especially  for  high  frequencies,  it  is  also  implied  that  near  regional 
range  body  waves  (out  200  km  or  so  from  a  source)  will  be  weakly  attenuated 
and  high  frequencies  should  be  propagated  efficiently  in  all  cases. 

In  regions  with  little  or  no  low  velocity  zones  (V  ) ,  the  range  of 
efficient  high  frequency  propagation  could  be  much  greater-perhaps 
out  to  15°  or  greater. 

Some  research  that  could  provide  needed  detail  and  better  quantify  the 
first  order  models  so  far  obtained,  includes: 

(1)  Simultaneous  matching  of  explosion  event  body  wave  seismograms 
in  the  near,  regional  and  teleseismic  distance  ranges  with  the  objective 
of  eliminating  uncertainties  in  the  source  function,  so  that  Q  models 
could  be  obtain  that  were  relatively  free  from  trade-off  problems  with  the 
source  function. 

/ 

(2)  Use  long  period  surface  waves  and  high  resolution  analysis 
methods  for  station  to  station  analysis  of  attenuation  to  obtain  Q 
models  that  would  be  free  from  source  trade-off  problems.  This  approach 
would  also  give  regional  Q  models. 
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Appendix  3 

Theoretical  vs  Mg  Curves  for  Earthquake  Sources 


/ 


A.  Theoretical  m^  vs  Mg  curves  for  Thrust,  Dip  Slip  and  Strike  slip  events 
in  Basin  and  Range  type  continental-teotonio  structure. 

Figure  1-a.  Thrust  (45°)  events  at  a  depth  of  10  km. 

Figure  1-b.  Thurst  (45°)  events  at  a  depth  of  15  km. 

Figure  1-c.  Thrust  (45°)  events  at  a  depth  of  25  km. 

Figure  1-d.  Thrust  (45°)  events  at  a  depth  of  35  km. 

Figure  1-e.  Thrust  (45°)  events  at  a  depth  of  45  km. 

Figure  2-a.  Dip  slip  (normal)  events  at  a  depth  of  10  km. 

Figure  2-b.  Dip  slip  (normal)  events  at  a  depth  of  15  km. 

Figure  2-c.  Dip  slip  (normal)  events  at  a  depth  of  25  km. 

Figure  2-d.  Dip  slip  (normal)  events  at  a  depth  of  35  km. 

Figure  2-e.  Dip  slip  (normal)  events  at  a  depth  of  45  km. 

Figure  3-a.  Strike-slip  events  at  a  depth  of  10  km. 

Figure  3-b.  Strike-slip  events  at  a  depth  of  15  km. 

Figure  3-c.  Strike-slip  events  at  a  depth  of  25  km. 

Figure  3-d.  Strike-slip  events  at  a  depth  of  35  km. 

Figure  3-e.  Strike-slip  events  at  a  depth  of  45  km. 
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B.  Theoretical  !»{,  Vs  Mg  curves  for  Thrust,  Dip  Slip  and  Strike  slip  events 
in  Island  Arc  type  oceanic  structure. 


Figure  4-a. 
Figure  4-b. 
Figure  4-c. 
Figure  4-d. 
Figure  4-e. 


Thrust  (45°)  events  at  a  depth  of  10  km. 
Thrust  (45°)  events  at  a  depth  of  15  km. 
Thrust  (45°)  events  at  a  depth  of  25  km. 
Thrust  (45°)  events  at  a  depth  of  35  km. 
Thrust  (45°)  events  at  a  depth  of  45  km. 


Figure  5-a.  Dip  slip  (normal)  events  at  a 
Figure  5-b.  Dip  slip  (normal)  events  at  a 
Figure  5-c.  Dip  slip  (normal)  events  at  a 
Figure  5-d.  Dip  slip  (normal)  events  at  a 
Figure  5-e.  Dip  slip  (normal)  events  at  a 
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Figure  4-d. 
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Figure  6-d. 


